1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n); where nis bigger 0
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
76 //-------------------------------------------------------
81 #include "Riostream.h"
82 #include <TClonesArray.h>
84 #include <TObjArray.h>
87 #include "AliComplexCluster.h"
88 #include "AliESDEvent.h"
89 #include "AliESDtrack.h"
90 #include "AliESDVertex.h"
94 #include "AliRunLoader.h"
95 #include "AliTPCClustersRow.h"
96 #include "AliTPCParam.h"
97 #include "AliTPCReconstructor.h"
98 #include "AliTPCpolyTrack.h"
99 #include "AliTPCreco.h"
100 #include "AliTPCseed.h"
101 #include "AliTPCtrackerMI.h"
102 #include "TStopwatch.h"
103 #include "AliTPCReconstructor.h"
105 #include "TTreeStream.h"
106 #include "AliAlignObj.h"
107 #include "AliTrackPointArray.h"
109 #include "AliTPCcalibDB.h"
110 #include "AliTPCTransform.h"
114 ClassImp(AliTPCtrackerMI)
117 class AliTPCFastMath {
120 static Double_t FastAsin(Double_t x);
122 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
125 Double_t AliTPCFastMath::fgFastAsin[20000];
126 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
128 AliTPCFastMath::AliTPCFastMath(){
130 // initialized lookup table;
131 for (Int_t i=0;i<10000;i++){
132 fgFastAsin[2*i] = TMath::ASin(i/10000.);
133 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
137 Double_t AliTPCFastMath::FastAsin(Double_t x){
139 // return asin using lookup table
141 Int_t index = int(x*10000);
142 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
145 Int_t index = int(x*10000);
146 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
148 //__________________________________________________________________
149 AliTPCtrackerMI::AliTPCtrackerMI()
171 // default constructor
174 //_____________________________________________________________________
178 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
180 //update track information using current cluster - track->fCurrentCluster
183 AliTPCclusterMI* c =track->GetCurrentCluster();
184 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
186 UInt_t i = track->GetCurrentClusterIndex1();
188 Int_t sec=(i&0xff000000)>>24;
189 //Int_t row = (i&0x00ff0000)>>16;
190 track->SetRow((i&0x00ff0000)>>16);
191 track->SetSector(sec);
192 // Int_t index = i&0xFFFF;
193 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
194 track->SetClusterIndex2(track->GetRow(), i);
195 //track->fFirstPoint = row;
196 //if ( track->fLastPoint<row) track->fLastPoint =row;
197 // if (track->fRow<0 || track->fRow>160) {
198 // printf("problem\n");
200 if (track->GetFirstPoint()>track->GetRow())
201 track->SetFirstPoint(track->GetRow());
202 if (track->GetLastPoint()<track->GetRow())
203 track->SetLastPoint(track->GetRow());
206 track->SetClusterPointer(track->GetRow(),c);
209 Double_t angle2 = track->GetSnp()*track->GetSnp();
211 //SET NEW Track Point
213 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
215 angle2 = TMath::Sqrt(angle2/(1-angle2));
216 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
218 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
219 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
220 point.SetErrY(sqrt(track->GetErrorY2()));
221 point.SetErrZ(sqrt(track->GetErrorZ2()));
223 point.SetX(track->GetX());
224 point.SetY(track->GetY());
225 point.SetZ(track->GetZ());
226 point.SetAngleY(angle2);
227 point.SetAngleZ(track->GetTgl());
228 if (point.IsShared()){
229 track->SetErrorY2(track->GetErrorY2()*4);
230 track->SetErrorZ2(track->GetErrorZ2()*4);
234 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
236 track->SetErrorY2(track->GetErrorY2()*1.3);
237 track->SetErrorY2(track->GetErrorY2()+0.01);
238 track->SetErrorZ2(track->GetErrorZ2()*1.3);
239 track->SetErrorZ2(track->GetErrorZ2()+0.005);
241 if (accept>0) return 0;
242 if (track->GetNumberOfClusters()%20==0){
243 // if (track->fHelixIn){
244 // TClonesArray & larr = *(track->fHelixIn);
245 // Int_t ihelix = larr.GetEntriesFast();
246 // new(larr[ihelix]) AliHelix(*track) ;
249 track->SetNoCluster(0);
250 return track->Update(c,chi2,i);
255 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor,
256 Float_t cory, Float_t corz)
259 // decide according desired precision to accept given
260 // cluster for tracking
261 Double_t sy2=ErrY2(seed,cluster)*cory;
262 Double_t sz2=ErrZ2(seed,cluster)*corz;
263 //sy2=ErrY2(seed,cluster)*cory;
264 //sz2=ErrZ2(seed,cluster)*cory;
266 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
267 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
269 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
270 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
271 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
272 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
274 Double_t rdistance2 = rdistancey2+rdistancez2;
277 if (rdistance2>16) return 3;
280 if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)
281 return 2; //suspisiouce - will be changed
283 if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)
284 // strict cut on overlaped cluster
285 return 2; //suspisiouce - will be changed
287 if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor )
288 && cluster->GetType()<0){
289 seed->SetNFoundable(seed->GetNFoundable()-1);
298 //_____________________________________________________________________________
299 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
301 fkNIS(par->GetNInnerSector()/2),
303 fkNOS(par->GetNOuterSector()/2),
320 //---------------------------------------------------------------------
321 // The main TPC tracker constructor
322 //---------------------------------------------------------------------
323 fInnerSec=new AliTPCSector[fkNIS];
324 fOuterSec=new AliTPCSector[fkNOS];
327 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
328 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
331 Int_t nrowlow = par->GetNRowLow();
332 Int_t nrowup = par->GetNRowUp();
335 for (Int_t i=0;i<nrowlow;i++){
336 fXRow[i] = par->GetPadRowRadiiLow(i);
337 fPadLength[i]= par->GetPadPitchLength(0,i);
338 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
342 for (Int_t i=0;i<nrowup;i++){
343 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
344 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
345 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
348 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
350 //________________________________________________________________________
351 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
372 //------------------------------------
373 // dummy copy constructor
374 //------------------------------------------------------------------
377 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
378 //------------------------------
380 //--------------------------------------------------------------
383 //_____________________________________________________________________________
384 AliTPCtrackerMI::~AliTPCtrackerMI() {
385 //------------------------------------------------------------------
386 // TPC tracker destructor
387 //------------------------------------------------------------------
394 if (fDebugStreamer) delete fDebugStreamer;
397 void AliTPCtrackerMI::SetIO()
401 fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName());
403 fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName());
405 AliTPCtrack *iotrack= new AliTPCtrack;
406 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
412 void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESDEvent * event)
428 AliTPCtrack *iotrack= new AliTPCtrack;
429 // iotrack->fHelixIn = new TClonesArray("AliHelix");
430 //iotrack->fHelixOut = new TClonesArray("AliHelix");
431 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
434 if (output && (fDebug&2)){
435 //write the full seed information if specified in debug mode
437 fSeedTree = new TTree("Seeds","Seeds");
438 AliTPCseed * vseed = new AliTPCseed;
440 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
441 arrtr->ExpandCreateFast(160);
442 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
444 vseed->SetPoints(arrtr);
445 vseed->SetEPoints(arre);
446 // vseed->fClusterPoints = arrcl;
447 fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
450 fTreeDebug = new TTree("trackDebug","trackDebug");
451 TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
452 fTreeDebug->Branch("debug",&arrd,32000,99);
460 void AliTPCtrackerMI::FillESD(TObjArray* arr)
464 //fill esds using updated tracks
466 // write tracks to the event
467 // store index of the track
468 Int_t nseed=arr->GetEntriesFast();
469 //FindKinks(arr,fEvent);
470 for (Int_t i=0; i<nseed; i++) {
471 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
474 // pt->PropagateTo(fParam->GetInnerRadiusLow());
475 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
476 pt->PropagateTo(fParam->GetInnerRadiusLow());
479 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
481 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
482 iotrack.SetTPCPoints(pt->GetPoints());
483 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
484 iotrack.SetV0Indexes(pt->GetV0Indexes());
485 // iotrack.SetTPCpid(pt->fTPCr);
486 //iotrack.SetTPCindex(i);
487 fEvent->AddTrack(&iotrack);
491 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
493 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
494 iotrack.SetTPCPoints(pt->GetPoints());
495 //iotrack.SetTPCindex(i);
496 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
497 iotrack.SetV0Indexes(pt->GetV0Indexes());
498 // iotrack.SetTPCpid(pt->fTPCr);
499 fEvent->AddTrack(&iotrack);
503 // short tracks - maybe decays
505 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
506 Int_t found,foundable,shared;
507 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
508 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
510 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
511 //iotrack.SetTPCindex(i);
512 iotrack.SetTPCPoints(pt->GetPoints());
513 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
514 iotrack.SetV0Indexes(pt->GetV0Indexes());
515 //iotrack.SetTPCpid(pt->fTPCr);
516 fEvent->AddTrack(&iotrack);
521 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
522 Int_t found,foundable,shared;
523 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
524 if (found<20) continue;
525 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
528 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
529 iotrack.SetTPCPoints(pt->GetPoints());
530 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
531 iotrack.SetV0Indexes(pt->GetV0Indexes());
532 //iotrack.SetTPCpid(pt->fTPCr);
533 //iotrack.SetTPCindex(i);
534 fEvent->AddTrack(&iotrack);
537 // short tracks - secondaties
539 if ( (pt->GetNumberOfClusters()>30) ) {
540 Int_t found,foundable,shared;
541 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
542 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
544 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
545 iotrack.SetTPCPoints(pt->GetPoints());
546 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
547 iotrack.SetV0Indexes(pt->GetV0Indexes());
548 //iotrack.SetTPCpid(pt->fTPCr);
549 //iotrack.SetTPCindex(i);
550 fEvent->AddTrack(&iotrack);
555 if ( (pt->GetNumberOfClusters()>15)) {
556 Int_t found,foundable,shared;
557 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
558 if (found<15) continue;
559 if (foundable<=0) continue;
560 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
561 if (float(found)/float(foundable)<0.8) continue;
564 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
565 iotrack.SetTPCPoints(pt->GetPoints());
566 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
567 iotrack.SetV0Indexes(pt->GetV0Indexes());
568 // iotrack.SetTPCpid(pt->fTPCr);
569 //iotrack.SetTPCindex(i);
570 fEvent->AddTrack(&iotrack);
575 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
578 void AliTPCtrackerMI::WriteTracks(TTree * tree)
581 // write tracks from seed array to selected tree
585 AliTPCtrack *iotrack= new AliTPCtrack;
586 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
591 void AliTPCtrackerMI::WriteTracks()
594 // write tracks to the given output tree -
595 // output specified with SetIO routine
602 AliTPCtrack *iotrack= 0;
603 Int_t nseed=fSeeds->GetEntriesFast();
604 //for (Int_t i=0; i<nseed; i++) {
605 // iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
606 // if (iotrack) break;
608 //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
609 TBranch * br = fOutput->GetBranch("tracks");
610 br->SetAddress(&iotrack);
612 for (Int_t i=0; i<nseed; i++) {
613 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
615 AliTPCtrack * track = new AliTPCtrack(*pt);
618 // br->SetAddress(&iotrack);
623 //fOutput->GetDirectory()->cd();
629 //write the full seed information if specified in debug mode
631 AliTPCseed * vseed = new AliTPCseed;
633 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
634 arrtr->ExpandCreateFast(160);
635 //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
636 //arrcl->ExpandCreateFast(160);
637 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
639 vseed->SetPoints(arrtr);
640 vseed->SetEPoints(arre);
641 // vseed->fClusterPoints = arrcl;
642 //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
643 TBranch * brseed = fSeedTree->GetBranch("seeds");
645 Int_t nseed=fSeeds->GetEntriesFast();
647 for (Int_t i=0; i<nseed; i++) {
648 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
650 pt->SetPoints(arrtr);
651 // pt->fClusterPoints = arrcl;
652 pt->SetEPoints(arre);
655 brseed->SetAddress(&vseed);
659 // pt->fClusterPoints = 0;
662 if (fTreeDebug) fTreeDebug->Write();
670 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
673 //seed->SetErrorY2(0.1);
675 //calculate look-up table at the beginning
676 static Bool_t ginit = kFALSE;
677 static Float_t gnoise1,gnoise2,gnoise3;
678 static Float_t ggg1[10000];
679 static Float_t ggg2[10000];
680 static Float_t ggg3[10000];
681 static Float_t glandau1[10000];
682 static Float_t glandau2[10000];
683 static Float_t glandau3[10000];
685 static Float_t gcor01[500];
686 static Float_t gcor02[500];
687 static Float_t gcorp[500];
692 for (Int_t i=1;i<500;i++){
693 Float_t rsigma = float(i)/100.;
694 gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
695 gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
696 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
700 for (Int_t i=3;i<10000;i++){
704 Float_t amp = float(i);
705 Float_t padlength =0.75;
706 gnoise1 = 0.0004/padlength;
707 Float_t nel = 0.268*amp;
708 Float_t nprim = 0.155*amp;
709 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
710 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
711 if (glandau1[i]>1) glandau1[i]=1;
712 glandau1[i]*=padlength*padlength/12.;
716 gnoise2 = 0.0004/padlength;
719 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
720 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
721 if (glandau2[i]>1) glandau2[i]=1;
722 glandau2[i]*=padlength*padlength/12.;
727 gnoise3 = 0.0004/padlength;
730 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
731 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
732 if (glandau3[i]>1) glandau3[i]=1;
733 glandau3[i]*=padlength*padlength/12.;
741 Int_t amp = int(TMath::Abs(cl->GetQ()));
743 seed->SetErrorY2(1.);
747 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
748 Int_t ctype = cl->GetType();
749 Float_t padlength= GetPadPitchLength(seed->GetRow());
750 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
751 angle2 = angle2/(1-angle2);
754 Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
757 if (fSectors==fInnerSec){
759 res = ggg1[amp]*z+glandau1[amp]*angle2;
760 if (ctype==0) res *= gcor01[rsigmay];
763 res*= gcorp[rsigmay];
769 res = ggg2[amp]*z+glandau2[amp]*angle2;
770 if (ctype==0) res *= gcor02[rsigmay];
773 res*= gcorp[rsigmay];
778 res = ggg3[amp]*z+glandau3[amp]*angle2;
779 if (ctype==0) res *= gcor02[rsigmay];
782 res*= gcorp[rsigmay];
789 res*=2.4; // overestimate error 2 times
796 seed->SetErrorY2(res);
804 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
807 //seed->SetErrorY2(0.1);
809 //calculate look-up table at the beginning
810 static Bool_t ginit = kFALSE;
811 static Float_t gnoise1,gnoise2,gnoise3;
812 static Float_t ggg1[10000];
813 static Float_t ggg2[10000];
814 static Float_t ggg3[10000];
815 static Float_t glandau1[10000];
816 static Float_t glandau2[10000];
817 static Float_t glandau3[10000];
819 static Float_t gcor01[1000];
820 static Float_t gcor02[1000];
821 static Float_t gcorp[1000];
826 for (Int_t i=1;i<1000;i++){
827 Float_t rsigma = float(i)/100.;
828 gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
829 gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
830 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
834 for (Int_t i=3;i<10000;i++){
838 Float_t amp = float(i);
839 Float_t padlength =0.75;
840 gnoise1 = 0.0004/padlength;
841 Float_t nel = 0.268*amp;
842 Float_t nprim = 0.155*amp;
843 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
844 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
845 if (glandau1[i]>1) glandau1[i]=1;
846 glandau1[i]*=padlength*padlength/12.;
850 gnoise2 = 0.0004/padlength;
853 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
854 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
855 if (glandau2[i]>1) glandau2[i]=1;
856 glandau2[i]*=padlength*padlength/12.;
861 gnoise3 = 0.0004/padlength;
864 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
865 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
866 if (glandau3[i]>1) glandau3[i]=1;
867 glandau3[i]*=padlength*padlength/12.;
875 Int_t amp = int(TMath::Abs(cl->GetQ()));
877 seed->SetErrorY2(1.);
881 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
882 Int_t ctype = cl->GetType();
883 Float_t padlength= GetPadPitchLength(seed->GetRow());
885 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
886 // if (angle2<0.6) angle2 = 0.6;
887 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
890 Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
893 if (fSectors==fInnerSec){
895 res = ggg1[amp]*z+glandau1[amp]*angle2;
896 if (ctype==0) res *= gcor01[rsigmaz];
899 res*= gcorp[rsigmaz];
905 res = ggg2[amp]*z+glandau2[amp]*angle2;
906 if (ctype==0) res *= gcor02[rsigmaz];
909 res*= gcorp[rsigmaz];
914 res = ggg3[amp]*z+glandau3[amp]*angle2;
915 if (ctype==0) res *= gcor02[rsigmaz];
918 res*= gcorp[rsigmaz];
927 if ((ctype<0) &&<70){
935 seed->SetErrorZ2(res);
942 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
945 //seed->SetErrorZ2(0.1);
949 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
951 Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
952 Int_t ctype = cl->GetType();
953 Float_t amp = TMath::Abs(cl->GetQ());
958 Float_t landau=2 ; //landau fluctuation part
959 Float_t gg=2; // gg fluctuation part
960 Float_t padlength= GetPadPitchLength(seed->GetX());
962 if (fSectors==fInnerSec){
963 snoise2 = 0.0004/padlength;
966 gg = (2+0.001*nel/(padlength*padlength))/nel;
967 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
968 if (landau>1) landau=1;
971 snoise2 = 0.0004/padlength;
974 gg = (2+0.0008*nel/(padlength*padlength))/nel;
975 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
976 if (landau>1) landau=1;
978 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
981 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
982 angle2 = TMath::Sqrt((1-angle2));
983 if (angle2<0.6) angle2 = 0.6;
986 Float_t angle = seed->GetTgl()/angle2;
987 Float_t angular = landau*angle*angle*padlength*padlength/12.;
988 Float_t res = sdiff + angular;
991 if ((ctype==0) && (fSectors ==fOuterSec))
992 res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
994 if ((ctype==0) && (fSectors ==fInnerSec))
995 res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
999 res*= TMath::Power(rsigmaz+0.5,1.5); //0.31+0.147*ctype;
1005 if ((ctype<0) &&<70){
1013 seed->SetErrorZ2(res);
1021 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
1023 //rotate to track "local coordinata
1024 Float_t x = seed->GetX();
1025 Float_t y = seed->GetY();
1026 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
1029 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
1030 if (!seed->Rotate(fSectors->GetAlpha()))
1032 } else if (y <-ymax) {
1033 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
1034 if (!seed->Rotate(-fSectors->GetAlpha()))
1042 //_____________________________________________________________________________
1043 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1044 Double_t x2,Double_t y2,
1045 Double_t x3,Double_t y3)
1047 //-----------------------------------------------------------------
1048 // Initial approximation of the track curvature
1049 //-----------------------------------------------------------------
1050 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1051 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1052 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1053 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1054 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1056 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1057 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1058 return -xr*yr/sqrt(xr*xr+yr*yr);
1063 //_____________________________________________________________________________
1064 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1065 Double_t x2,Double_t y2,
1066 Double_t x3,Double_t y3)
1068 //-----------------------------------------------------------------
1069 // Initial approximation of the track curvature
1070 //-----------------------------------------------------------------
1076 Double_t det = x3*y2-x2*y3;
1081 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1082 Double_t x0 = x3*0.5-y3*u;
1083 Double_t y0 = y3*0.5+x3*u;
1084 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1090 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1091 Double_t x2,Double_t y2,
1092 Double_t x3,Double_t y3)
1094 //-----------------------------------------------------------------
1095 // Initial approximation of the track curvature
1096 //-----------------------------------------------------------------
1102 Double_t det = x3*y2-x2*y3;
1107 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1108 Double_t x0 = x3*0.5-y3*u;
1109 Double_t y0 = y3*0.5+x3*u;
1110 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1119 //_____________________________________________________________________________
1120 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1121 Double_t x2,Double_t y2,
1122 Double_t x3,Double_t y3)
1124 //-----------------------------------------------------------------
1125 // Initial approximation of the track curvature times center of curvature
1126 //-----------------------------------------------------------------
1127 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1128 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1129 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1130 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1131 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1133 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1135 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1138 //_____________________________________________________________________________
1139 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1140 Double_t x2,Double_t y2,
1141 Double_t z1,Double_t z2)
1143 //-----------------------------------------------------------------
1144 // Initial approximation of the tangent of the track dip angle
1145 //-----------------------------------------------------------------
1146 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1150 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1151 Double_t x2,Double_t y2,
1152 Double_t z1,Double_t z2, Double_t c)
1154 //-----------------------------------------------------------------
1155 // Initial approximation of the tangent of the track dip angle
1156 //-----------------------------------------------------------------
1160 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1162 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1163 if (TMath::Abs(d*c*0.5)>1) return 0;
1164 // Double_t angle2 = TMath::ASin(d*c*0.5);
1165 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1166 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1168 angle2 = (z1-z2)*c/(angle2*2.);
1172 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1173 {//-----------------------------------------------------------------
1174 // This function find proloncation of a track to a reference plane x=x2.
1175 //-----------------------------------------------------------------
1179 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1183 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1184 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1188 Double_t dy = dx*(c1+c2)/(r1+r2);
1191 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1193 if (TMath::Abs(delta)>0.01){
1194 dz = x[3]*TMath::ASin(delta)/x[4];
1196 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1199 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1207 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1212 return LoadClusters();
1215 Int_t AliTPCtrackerMI::LoadClusters()
1218 // load clusters to the memory
1219 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1220 clrow->SetClass("AliTPCclusterMI");
1222 clrow->GetArray()->ExpandCreateFast(10000);
1224 // TTree * tree = fClustersArray.GetTree();
1226 TTree * tree = fInput;
1227 TBranch * br = tree->GetBranch("Segment");
1228 br->SetAddress(&clrow);
1230 Int_t j=Int_t(tree->GetEntries());
1231 for (Int_t i=0; i<j; i++) {
1235 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1236 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1237 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1240 AliTPCRow * tpcrow=0;
1243 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1247 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1248 left = (sec-fkNIS*2)/fkNOS;
1251 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1252 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1253 for (Int_t i=0;i<tpcrow->GetN1();i++)
1254 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1257 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1258 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1259 for (Int_t i=0;i<tpcrow->GetN2();i++)
1260 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1271 void AliTPCtrackerMI::UnloadClusters()
1274 // unload clusters from the memory
1276 Int_t nrows = fOuterSec->GetNRows();
1277 for (Int_t sec = 0;sec<fkNOS;sec++)
1278 for (Int_t row = 0;row<nrows;row++){
1279 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1281 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1282 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1284 tpcrow->ResetClusters();
1287 nrows = fInnerSec->GetNRows();
1288 for (Int_t sec = 0;sec<fkNIS;sec++)
1289 for (Int_t row = 0;row<nrows;row++){
1290 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1292 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1293 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1295 tpcrow->ResetClusters();
1301 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1305 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1307 AliFatal("Tranformations not in calibDB");
1309 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1310 Int_t i[1]={cluster->GetDetector()};
1311 transform->Transform(x,i,0,1);
1312 if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
1313 if (cluster->GetDetector()%36>17){
1319 // in debug mode check the transformation
1321 if (AliTPCReconstructor::StreamLevel()>1) {
1322 TTreeSRedirector &cstream = *fDebugStreamer;
1323 cstream<<"Transform"<<
1330 cluster->SetX(x[0]);
1331 cluster->SetY(x[1]);
1332 cluster->SetZ(x[2]);
1338 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1339 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1341 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1342 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1343 if (mat) mat->LocalToMaster(pos,posC);
1345 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1347 cluster->SetX(posC[0]);
1348 cluster->SetY(posC[1]);
1349 cluster->SetZ(posC[2]);
1352 //_____________________________________________________________________________
1353 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1354 //-----------------------------------------------------------------
1355 // This function fills outer TPC sectors with clusters.
1356 //-----------------------------------------------------------------
1357 Int_t nrows = fOuterSec->GetNRows();
1359 for (Int_t sec = 0;sec<fkNOS;sec++)
1360 for (Int_t row = 0;row<nrows;row++){
1361 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1362 Int_t sec2 = sec+2*fkNIS;
1364 Int_t ncl = tpcrow->GetN1();
1366 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1367 index=(((sec2<<8)+row)<<16)+ncl;
1368 tpcrow->InsertCluster(c,index);
1371 ncl = tpcrow->GetN2();
1373 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1374 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1375 tpcrow->InsertCluster(c,index);
1378 // write indexes for fast acces
1380 for (Int_t i=0;i<510;i++)
1381 tpcrow->SetFastCluster(i,-1);
1382 for (Int_t i=0;i<tpcrow->GetN();i++){
1383 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1384 tpcrow->SetFastCluster(zi,i); // write index
1387 for (Int_t i=0;i<510;i++){
1388 if (tpcrow->GetFastCluster(i)<0)
1389 tpcrow->SetFastCluster(i,last);
1391 last = tpcrow->GetFastCluster(i);
1400 //_____________________________________________________________________________
1401 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1402 //-----------------------------------------------------------------
1403 // This function fills inner TPC sectors with clusters.
1404 //-----------------------------------------------------------------
1405 Int_t nrows = fInnerSec->GetNRows();
1407 for (Int_t sec = 0;sec<fkNIS;sec++)
1408 for (Int_t row = 0;row<nrows;row++){
1409 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1412 Int_t ncl = tpcrow->GetN1();
1414 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1415 index=(((sec<<8)+row)<<16)+ncl;
1416 tpcrow->InsertCluster(c,index);
1419 ncl = tpcrow->GetN2();
1421 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1422 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1423 tpcrow->InsertCluster(c,index);
1426 // write indexes for fast acces
1428 for (Int_t i=0;i<510;i++)
1429 tpcrow->SetFastCluster(i,-1);
1430 for (Int_t i=0;i<tpcrow->GetN();i++){
1431 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1432 tpcrow->SetFastCluster(zi,i); // write index
1435 for (Int_t i=0;i<510;i++){
1436 if (tpcrow->GetFastCluster(i)<0)
1437 tpcrow->SetFastCluster(i,last);
1439 last = tpcrow->GetFastCluster(i);
1451 //_________________________________________________________________________
1452 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1453 //--------------------------------------------------------------------
1454 // Return pointer to a given cluster
1455 //--------------------------------------------------------------------
1456 if (index<0) return 0; // no cluster
1457 Int_t sec=(index&0xff000000)>>24;
1458 Int_t row=(index&0x00ff0000)>>16;
1459 Int_t ncl=(index&0x00007fff)>>00;
1461 const AliTPCRow * tpcrow=0;
1462 AliTPCclusterMI * clrow =0;
1464 if (sec<0 || sec>=fkNIS*4) {
1465 AliWarning(Form("Wrong sector %d",sec));
1470 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1471 if (tpcrow==0) return 0;
1474 if (tpcrow->GetN1()<=ncl) return 0;
1475 clrow = tpcrow->GetClusters1();
1478 if (tpcrow->GetN2()<=ncl) return 0;
1479 clrow = tpcrow->GetClusters2();
1483 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1484 if (tpcrow==0) return 0;
1486 if (sec-2*fkNIS<fkNOS) {
1487 if (tpcrow->GetN1()<=ncl) return 0;
1488 clrow = tpcrow->GetClusters1();
1491 if (tpcrow->GetN2()<=ncl) return 0;
1492 clrow = tpcrow->GetClusters2();
1496 return &(clrow[ncl]);
1502 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1503 //-----------------------------------------------------------------
1504 // This function tries to find a track prolongation to next pad row
1505 //-----------------------------------------------------------------
1507 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1508 AliTPCclusterMI *cl=0;
1509 Int_t tpcindex= t.GetClusterIndex2(nr);
1511 // update current shape info every 5 pad-row
1512 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1516 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1518 if (tpcindex==-1) return 0; //track in dead zone
1520 cl = t.GetClusterPointer(nr);
1521 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1522 t.SetCurrentClusterIndex1(tpcindex);
1525 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1526 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1528 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1529 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1531 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1532 Double_t rotation = angle-t.GetAlpha();
1533 t.SetRelativeSector(relativesector);
1534 if (!t.Rotate(rotation)) return 0;
1536 if (!t.PropagateTo(x)) return 0;
1538 t.SetCurrentCluster(cl);
1540 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1541 if ((tpcindex&0x8000)==0) accept =0;
1543 //if founded cluster is acceptible
1544 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1545 t.SetErrorY2(t.GetErrorY2()+0.03);
1546 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1547 t.SetErrorY2(t.GetErrorY2()*3);
1548 t.SetErrorZ2(t.GetErrorZ2()*3);
1550 t.SetNFoundable(t.GetNFoundable()+1);
1551 UpdateTrack(&t,accept);
1556 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1558 // not look for new cluster during refitting
1559 t.SetNFoundable(t.GetNFoundable()+1);
1564 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1565 Double_t y=t.GetYat(x);
1566 if (TMath::Abs(y)>ymax){
1568 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1569 if (!t.Rotate(fSectors->GetAlpha()))
1571 } else if (y <-ymax) {
1572 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1573 if (!t.Rotate(-fSectors->GetAlpha()))
1579 if (!t.PropagateTo(x)) {
1580 if (fIteration==0) t.SetRemoval(10);
1584 Double_t z=t.GetZ();
1587 if (!IsActive(t.GetRelativeSector(),nr)) {
1589 t.SetClusterIndex2(nr,-1);
1592 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1593 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1594 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1596 if (!isActive || !isActive2) return 0;
1598 const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1599 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1601 Double_t roadz = 1.;
1603 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1605 t.SetClusterIndex2(nr,-1);
1610 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1611 t.SetNFoundable(t.GetNFoundable()+1);
1617 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1618 cl = krow.FindNearest2(y,z,roady,roadz,index);
1619 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1622 t.SetCurrentCluster(cl);
1624 if (fIteration==2&&cl->IsUsed(10)) return 0;
1625 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1626 if (fIteration==2&&cl->IsUsed(11)) {
1627 t.SetErrorY2(t.GetErrorY2()+0.03);
1628 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1629 t.SetErrorY2(t.GetErrorY2()*3);
1630 t.SetErrorZ2(t.GetErrorZ2()*3);
1633 if (t.fCurrentCluster->IsUsed(10)){
1638 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1644 if (accept<3) UpdateTrack(&t,accept);
1647 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1653 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1654 //-----------------------------------------------------------------
1655 // This function tries to find a track prolongation to next pad row
1656 //-----------------------------------------------------------------
1658 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1660 if (!t.GetProlongation(x,y,z)) {
1666 if (TMath::Abs(y)>ymax){
1669 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1670 if (!t.Rotate(fSectors->GetAlpha()))
1672 } else if (y <-ymax) {
1673 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1674 if (!t.Rotate(-fSectors->GetAlpha()))
1677 if (!t.PropagateTo(x)) {
1680 t.GetProlongation(x,y,z);
1683 // update current shape info every 3 pad-row
1684 if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1685 // t.fCurrentSigmaY = GetSigmaY(&t);
1686 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1690 AliTPCclusterMI *cl=0;
1695 const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1696 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1698 Double_t roadz = 1.;
1701 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1703 t.SetClusterIndex2(row,-1);
1708 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1712 if ((cl==0)&&(krow)) {
1713 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1714 cl = krow.FindNearest2(y,z,roady,roadz,index);
1716 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1720 t.SetCurrentCluster(cl);
1721 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1723 t.SetClusterIndex2(row,index);
1724 t.SetClusterPointer(row, cl);
1732 //_________________________________________________________________________
1733 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1735 // Get track space point by index
1736 // return false in case the cluster doesn't exist
1737 AliTPCclusterMI *cl = GetClusterMI(index);
1738 if (!cl) return kFALSE;
1739 Int_t sector = (index&0xff000000)>>24;
1740 // Int_t row = (index&0x00ff0000)>>16;
1742 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1743 xyz[0] = cl->GetX();
1744 xyz[1] = cl->GetY();
1745 xyz[2] = cl->GetZ();
1747 fParam->AdjustCosSin(sector,cos,sin);
1748 Float_t x = cos*xyz[0]-sin*xyz[1];
1749 Float_t y = cos*xyz[1]+sin*xyz[0];
1751 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1752 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1753 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1754 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1755 cov[0] = sin*sin*sigmaY2;
1756 cov[1] = -sin*cos*sigmaY2;
1758 cov[3] = cos*cos*sigmaY2;
1761 p.SetXYZ(x,y,xyz[2],cov);
1762 AliGeomManager::ELayerID iLayer;
1764 if (sector < fParam->GetNInnerSector()) {
1765 iLayer = AliGeomManager::kTPC1;
1769 iLayer = AliGeomManager::kTPC2;
1770 idet = sector - fParam->GetNInnerSector();
1772 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1773 p.SetVolumeID(volid);
1779 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1780 //-----------------------------------------------------------------
1781 // This function tries to find a track prolongation to next pad row
1782 //-----------------------------------------------------------------
1783 t.SetCurrentCluster(0);
1784 t.SetCurrentClusterIndex1(0);
1786 Double_t xt=t.GetX();
1787 Int_t row = GetRowNumber(xt)-1;
1788 Double_t ymax= GetMaxY(nr);
1790 if (row < nr) return 1; // don't prolongate if not information until now -
1791 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1793 // return 0; // not prolongate strongly inclined tracks
1795 // if (TMath::Abs(t.GetSnp())>0.95) {
1797 // return 0; // not prolongate strongly inclined tracks
1798 // }// patch 28 fev 06
1800 Double_t x= GetXrow(nr);
1802 //t.PropagateTo(x+0.02);
1803 //t.PropagateTo(x+0.01);
1804 if (!t.PropagateTo(x)){
1811 if (TMath::Abs(y)>ymax){
1813 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1814 if (!t.Rotate(fSectors->GetAlpha()))
1816 } else if (y <-ymax) {
1817 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1818 if (!t.Rotate(-fSectors->GetAlpha()))
1821 // if (!t.PropagateTo(x)){
1828 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1830 if (!IsActive(t.GetRelativeSector(),nr)) {
1832 t.SetClusterIndex2(nr,-1);
1835 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1837 AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1839 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1841 t.SetClusterIndex2(nr,-1);
1846 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1852 if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1853 // t.fCurrentSigmaY = GetSigmaY(&t);
1854 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1858 AliTPCclusterMI *cl=0;
1861 Double_t roady = 1.;
1862 Double_t roadz = 1.;
1866 index = t.GetClusterIndex2(nr);
1867 if ( (index>0) && (index&0x8000)==0){
1868 cl = t.GetClusterPointer(nr);
1869 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1870 t.SetCurrentClusterIndex1(index);
1872 t.SetCurrentCluster(cl);
1878 // if (index<0) return 0;
1879 UInt_t uindex = TMath::Abs(index);
1882 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1883 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1886 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1887 t.SetCurrentCluster(cl);
1893 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1894 //-----------------------------------------------------------------
1895 // This function tries to find a track prolongation to next pad row
1896 //-----------------------------------------------------------------
1898 //update error according neighborhoud
1900 if (t.GetCurrentCluster()) {
1902 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1904 if (t.GetCurrentCluster()->IsUsed(10)){
1909 t.SetNShared(t.GetNShared()+1);
1910 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1915 if (fIteration>0) accept = 0;
1916 if (accept<3) UpdateTrack(&t,accept);
1920 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1921 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1923 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1931 //_____________________________________________________________________________
1932 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1933 //-----------------------------------------------------------------
1934 // This function tries to find a track prolongation.
1935 //-----------------------------------------------------------------
1936 Double_t xt=t.GetX();
1938 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1939 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1940 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1942 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1944 Int_t first = GetRowNumber(xt)-1;
1945 for (Int_t nr= first; nr>=rf; nr-=step) {
1947 if (t.GetKinkIndexes()[0]>0){
1948 for (Int_t i=0;i<3;i++){
1949 Int_t index = t.GetKinkIndexes()[i];
1950 if (index==0) break;
1951 if (index<0) continue;
1953 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1955 printf("PROBLEM\n");
1958 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1960 AliExternalTrackParam paramd(t);
1961 kink->SetDaughter(paramd);
1962 kink->SetStatus(2,5);
1969 if (nr==80) t.UpdateReference();
1970 if (nr<fInnerSec->GetNRows())
1971 fSectors = fInnerSec;
1973 fSectors = fOuterSec;
1974 if (FollowToNext(t,nr)==0)
1983 //_____________________________________________________________________________
1984 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1985 //-----------------------------------------------------------------
1986 // This function tries to find a track prolongation.
1987 //-----------------------------------------------------------------
1988 Double_t xt=t.GetX();
1990 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1991 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1992 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1993 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1995 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1997 if (FollowToNextFast(t,nr)==0)
1998 if (!t.IsActive()) return 0;
2008 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
2009 //-----------------------------------------------------------------
2010 // This function tries to find a track prolongation.
2011 //-----------------------------------------------------------------
2013 Double_t xt=t.GetX();
2014 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
2015 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2016 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2017 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2019 Int_t first = t.GetFirstPoint();
2020 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
2022 if (first<0) first=0;
2023 for (Int_t nr=first; nr<=rf; nr++) {
2024 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2025 if (t.GetKinkIndexes()[0]<0){
2026 for (Int_t i=0;i<3;i++){
2027 Int_t index = t.GetKinkIndexes()[i];
2028 if (index==0) break;
2029 if (index>0) continue;
2030 index = TMath::Abs(index);
2031 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2033 printf("PROBLEM\n");
2036 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2038 AliExternalTrackParam paramm(t);
2039 kink->SetMother(paramm);
2040 kink->SetStatus(2,1);
2047 if (nr<fInnerSec->GetNRows())
2048 fSectors = fInnerSec;
2050 fSectors = fOuterSec;
2061 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2069 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2072 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2074 Float_t distance = TMath::Sqrt(dz2+dy2);
2075 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2078 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2079 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2084 if (firstpoint>lastpoint) {
2085 firstpoint =lastpoint;
2090 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2091 if (s1->GetClusterIndex2(i)>0) sum1++;
2092 if (s2->GetClusterIndex2(i)>0) sum2++;
2093 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2097 if (sum<5) return 0;
2099 Float_t summin = TMath::Min(sum1+1,sum2+1);
2100 Float_t ratio = (sum+1)/Float_t(summin);
2104 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2108 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2109 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2110 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2111 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2116 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2117 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2118 Int_t firstpoint = 0;
2119 Int_t lastpoint = 160;
2121 // if (firstpoint>=lastpoint-5) return;;
2123 for (Int_t i=firstpoint;i<lastpoint;i++){
2124 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2125 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2127 s1->SetSharedMapBit(i, kTRUE);
2128 s2->SetSharedMapBit(i, kTRUE);
2130 if (s1->GetClusterIndex2(i)>0)
2131 s1->SetClusterMapBit(i, kTRUE);
2133 if (sumshared>cutN0){
2136 for (Int_t i=firstpoint;i<lastpoint;i++){
2137 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2138 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2139 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2140 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2141 if (s1->IsActive()&&s2->IsActive()){
2142 p1->SetShared(kTRUE);
2143 p2->SetShared(kTRUE);
2149 if (sumshared>cutN0){
2150 for (Int_t i=0;i<4;i++){
2151 if (s1->GetOverlapLabel(3*i)==0){
2152 s1->SetOverlapLabel(3*i, s2->GetLabel());
2153 s1->SetOverlapLabel(3*i+1,sumshared);
2154 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2158 for (Int_t i=0;i<4;i++){
2159 if (s2->GetOverlapLabel(3*i)==0){
2160 s2->SetOverlapLabel(3*i, s1->GetLabel());
2161 s2->SetOverlapLabel(3*i+1,sumshared);
2162 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2169 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2172 //sort trackss according sectors
2174 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2175 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2177 //if (pt) RotateToLocal(pt);
2181 arr->Sort(); // sorting according relative sectors
2182 arr->Expand(arr->GetEntries());
2185 Int_t nseed=arr->GetEntriesFast();
2186 for (Int_t i=0; i<nseed; i++) {
2187 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2189 for (Int_t j=0;j<=12;j++){
2190 pt->SetOverlapLabel(j,0);
2193 for (Int_t i=0; i<nseed; i++) {
2194 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2196 if (pt->GetRemoval()>10) continue;
2197 for (Int_t j=i+1; j<nseed; j++){
2198 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2199 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2201 if (pt2->GetRemoval()<=10) {
2202 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2210 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2213 //sort tracks in array according mode criteria
2214 Int_t nseed = arr->GetEntriesFast();
2215 for (Int_t i=0; i<nseed; i++) {
2216 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2227 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2230 // Loop over all tracks and remove overlaped tracks (with lower quality)
2232 // 1. Unsign clusters
2233 // 2. Sort tracks according quality
2234 // Quality is defined by the number of cluster between first and last points
2236 // 3. Loop over tracks - decreasing quality order
2237 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2238 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2239 // c.) if track accepted - sign clusters
2241 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2242 // - AliTPCtrackerMI::PropagateBack()
2243 // - AliTPCtrackerMI::RefitInward()
2249 Int_t nseed = arr->GetEntriesFast();
2250 Float_t * quality = new Float_t[nseed];
2251 Int_t * indexes = new Int_t[nseed];
2255 for (Int_t i=0; i<nseed; i++) {
2256 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2261 pt->UpdatePoints(); //select first last max dens points
2262 Float_t * points = pt->GetPoints();
2263 if (points[3]<0.8) quality[i] =-1;
2264 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2265 //prefer high momenta tracks if overlaps
2266 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2268 TMath::Sort(nseed,quality,indexes);
2271 for (Int_t itrack=0; itrack<nseed; itrack++) {
2272 Int_t trackindex = indexes[itrack];
2273 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2276 if (quality[trackindex]<0){
2278 delete arr->RemoveAt(trackindex);
2281 arr->RemoveAt(trackindex);
2287 Int_t first = Int_t(pt->GetPoints()[0]);
2288 Int_t last = Int_t(pt->GetPoints()[2]);
2289 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2291 Int_t found,foundable,shared;
2292 pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE); // better to get statistic in "high-dens" region do't use full track as in line bellow
2293 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2294 Bool_t itsgold =kFALSE;
2297 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2301 if (Float_t(shared+1)/Float_t(found+1)>factor){
2302 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2303 delete arr->RemoveAt(trackindex);
2306 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2307 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2308 delete arr->RemoveAt(trackindex);
2314 //if (sharedfactor>0.4) continue;
2315 if (pt->GetKinkIndexes()[0]>0) continue;
2316 //Remove tracks with undefined properties - seems
2317 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2319 for (Int_t i=first; i<last; i++) {
2320 Int_t index=pt->GetClusterIndex2(i);
2321 // if (index<0 || index&0x8000 ) continue;
2322 if (index<0 || index&0x8000 ) continue;
2323 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2330 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2336 void AliTPCtrackerMI::UnsignClusters()
2339 // loop over all clusters and unsign them
2342 for (Int_t sec=0;sec<fkNIS;sec++){
2343 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2344 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2345 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2346 // if (cl[icl].IsUsed(10))
2348 cl = fInnerSec[sec][row].GetClusters2();
2349 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2350 //if (cl[icl].IsUsed(10))
2355 for (Int_t sec=0;sec<fkNOS;sec++){
2356 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2357 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2358 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2359 //if (cl[icl].IsUsed(10))
2361 cl = fOuterSec[sec][row].GetClusters2();
2362 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2363 //if (cl[icl].IsUsed(10))
2372 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2375 //sign clusters to be "used"
2377 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2378 // loop over "primaries"
2392 Int_t nseed = arr->GetEntriesFast();
2393 for (Int_t i=0; i<nseed; i++) {
2394 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2398 if (!(pt->IsActive())) continue;
2399 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2400 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2402 sumdens2+= dens*dens;
2403 sumn += pt->GetNumberOfClusters();
2404 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2405 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2408 sumchi2 +=chi2*chi2;
2413 Float_t mdensity = 0.9;
2414 Float_t meann = 130;
2415 Float_t meanchi = 1;
2416 Float_t sdensity = 0.1;
2417 Float_t smeann = 10;
2418 Float_t smeanchi =0.4;
2422 mdensity = sumdens/sum;
2424 meanchi = sumchi/sum;
2426 sdensity = sumdens2/sum-mdensity*mdensity;
2428 sdensity = TMath::Sqrt(sdensity);
2432 smeann = sumn2/sum-meann*meann;
2434 smeann = TMath::Sqrt(smeann);
2438 smeanchi = sumchi2/sum - meanchi*meanchi;
2440 smeanchi = TMath::Sqrt(smeanchi);
2446 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2448 for (Int_t i=0; i<nseed; i++) {
2449 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2453 if (pt->GetBSigned()) continue;
2454 if (pt->GetBConstrain()) continue;
2455 //if (!(pt->IsActive())) continue;
2457 Int_t found,foundable,shared;
2458 pt->GetClusterStatistic(0,160,found, foundable,shared);
2459 if (shared/float(found)>0.3) {
2460 if (shared/float(found)>0.9 ){
2461 //delete arr->RemoveAt(i);
2466 Bool_t isok =kFALSE;
2467 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2469 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2471 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2473 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2477 for (Int_t i=0; i<160; i++) {
2478 Int_t index=pt->GetClusterIndex2(i);
2479 if (index<0) continue;
2480 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2482 //if (!(c->IsUsed(10))) c->Use();
2489 Double_t maxchi = meanchi+2.*smeanchi;
2491 for (Int_t i=0; i<nseed; i++) {
2492 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2496 //if (!(pt->IsActive())) continue;
2497 if (pt->GetBSigned()) continue;
2498 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2499 if (chi>maxchi) continue;
2502 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2504 //sign only tracks with enoug big density at the beginning
2506 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2509 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2510 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2512 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2513 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2516 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2517 //Int_t noc=pt->GetNumberOfClusters();
2518 pt->SetBSigned(kTRUE);
2519 for (Int_t i=0; i<160; i++) {
2521 Int_t index=pt->GetClusterIndex2(i);
2522 if (index<0) continue;
2523 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2525 // if (!(c->IsUsed(10))) c->Use();
2530 // gLastCheck = nseed;
2538 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2540 // stop not active tracks
2541 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2542 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2543 Int_t nseed = arr->GetEntriesFast();
2545 for (Int_t i=0; i<nseed; i++) {
2546 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2550 if (!(pt->IsActive())) continue;
2551 StopNotActive(pt,row0,th0, th1,th2);
2557 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2560 // stop not active tracks
2561 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2562 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2565 Int_t foundable = 0;
2566 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2567 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2568 seed->Desactivate(10) ;
2572 for (Int_t i=row0; i<maxindex; i++){
2573 Int_t index = seed->GetClusterIndex2(i);
2574 if (index!=-1) foundable++;
2576 if (foundable<=30) sumgood1++;
2577 if (foundable<=50) {
2584 if (foundable>=30.){
2585 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2588 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2592 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2595 // back propagation of ESD tracks
2598 const Int_t kMaxFriendTracks=2000;
2602 //PrepareForProlongation(fSeeds,1);
2603 PropagateForward2(fSeeds);
2604 RemoveUsed2(fSeeds,0.4,0.4,20);
2606 TObjArray arraySeed(fSeeds->GetEntries());
2607 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2608 arraySeed.AddAt(fSeeds->At(i),i);
2610 SignShared(&arraySeed);
2611 // FindCurling(fSeeds, event,2); // find multi found tracks
2612 FindSplitted(fSeeds, event,2); // find multi found tracks
2615 Int_t nseed = fSeeds->GetEntriesFast();
2616 for (Int_t i=0;i<nseed;i++){
2617 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2618 if (!seed) continue;
2619 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2621 seed->PropagateTo(fParam->GetInnerRadiusLow());
2622 seed->UpdatePoints();
2624 AliESDtrack *esd=event->GetTrack(i);
2625 seed->CookdEdx(0.02,0.6);
2626 CookLabel(seed,0.1); //For comparison only
2627 esd->SetTPCClusterMap(seed->GetClusterMap());
2628 esd->SetTPCSharedMap(seed->GetSharedMap());
2630 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2631 TTreeSRedirector &cstream = *fDebugStreamer;
2638 if (seed->GetNumberOfClusters()>15){
2639 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2640 esd->SetTPCPoints(seed->GetPoints());
2641 esd->SetTPCPointsF(seed->GetNFoundable());
2642 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2643 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2644 Float_t dedx = seed->GetdEdx();
2645 esd->SetTPCsignal(dedx, sdedx, ndedx);
2647 // add seed to the esd track in Calib level
2649 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2650 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2651 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2652 esd->AddCalibObject(seedCopy);
2657 //printf("problem\n");
2660 //FindKinks(fSeeds,event);
2661 Info("RefitInward","Number of refitted tracks %d",ntracks);
2668 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2671 // back propagation of ESD tracks
2677 PropagateBack(fSeeds);
2678 RemoveUsed2(fSeeds,0.4,0.4,20);
2679 //FindCurling(fSeeds, fEvent,1);
2680 FindSplitted(fSeeds, event,1); // find multi found tracks
2683 Int_t nseed = fSeeds->GetEntriesFast();
2685 for (Int_t i=0;i<nseed;i++){
2686 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2687 if (!seed) continue;
2688 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2689 seed->UpdatePoints();
2690 AliESDtrack *esd=event->GetTrack(i);
2691 seed->CookdEdx(0.02,0.6);
2692 CookLabel(seed,0.1); //For comparison only
2693 if (seed->GetNumberOfClusters()>15){
2694 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2695 esd->SetTPCPoints(seed->GetPoints());
2696 esd->SetTPCPointsF(seed->GetNFoundable());
2697 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2698 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2699 Float_t dedx = seed->GetdEdx();
2700 esd->SetTPCsignal(dedx, sdedx, ndedx);
2702 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2703 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2704 if (AliTPCReconstructor::StreamLevel()>1) {
2705 (*fDebugStreamer)<<"Cback"<<
2707 "EventNrInFile="<<eventnumber<<
2708 "\n"; // patch 28 fev 06
2712 //FindKinks(fSeeds,event);
2713 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2721 void AliTPCtrackerMI::DeleteSeeds()
2726 Int_t nseed = fSeeds->GetEntriesFast();
2727 for (Int_t i=0;i<nseed;i++){
2728 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2729 if (seed) delete fSeeds->RemoveAt(i);
2736 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2739 //read seeds from the event
2741 Int_t nentr=event->GetNumberOfTracks();
2743 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2748 fSeeds = new TObjArray(nentr);
2752 for (Int_t i=0; i<nentr; i++) {
2753 AliESDtrack *esd=event->GetTrack(i);
2754 ULong_t status=esd->GetStatus();
2755 if (!(status&AliESDtrack::kTPCin)) continue;
2756 AliTPCtrack t(*esd);
2757 t.SetNumberOfClusters(0);
2758 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2759 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2760 seed->SetUniqueID(esd->GetID());
2761 for (Int_t ikink=0;ikink<3;ikink++) {
2762 Int_t index = esd->GetKinkIndex(ikink);
2763 seed->GetKinkIndexes()[ikink] = index;
2764 if (index==0) continue;
2765 index = TMath::Abs(index);
2766 AliESDkink * kink = fEvent->GetKink(index-1);
2767 if (kink&&esd->GetKinkIndex(ikink)<0){
2768 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2769 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2771 if (kink&&esd->GetKinkIndex(ikink)>0){
2772 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2773 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2777 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2778 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2779 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2784 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2785 Double_t par0[5],par1[5],alpha,x;
2786 esd->GetInnerExternalParameters(alpha,x,par0);
2787 esd->GetExternalParameters(x,par1);
2788 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2789 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2791 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2792 //reset covariance if suspicious
2793 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2794 seed->ResetCovariance(10.);
2799 // rotate to the local coordinate system
2801 fSectors=fInnerSec; fN=fkNIS;
2802 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2803 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2804 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2805 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2806 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2807 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2808 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2809 alpha-=seed->GetAlpha();
2810 if (!seed->Rotate(alpha)) {
2816 if (esd->GetKinkIndex(0)<=0){
2817 for (Int_t irow=0;irow<160;irow++){
2818 Int_t index = seed->GetClusterIndex2(irow);
2821 AliTPCclusterMI * cl = GetClusterMI(index);
2822 seed->SetClusterPointer(irow,cl);
2824 if ((index & 0x8000)==0){
2825 cl->Use(10); // accepted cluster
2827 cl->Use(6); // close cluster not accepted
2830 Info("ReadSeeds","Not found cluster");
2835 fSeeds->AddAt(seed,i);
2841 //_____________________________________________________________________________
2842 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2843 Float_t deltay, Int_t ddsec) {
2844 //-----------------------------------------------------------------
2845 // This function creates track seeds.
2846 // SEEDING WITH VERTEX CONSTRAIN
2847 //-----------------------------------------------------------------
2848 // cuts[0] - fP4 cut
2849 // cuts[1] - tan(phi) cut
2850 // cuts[2] - zvertex cut
2851 // cuts[3] - fP3 cut
2859 Double_t x[5], c[15];
2860 // Int_t di = i1-i2;
2862 AliTPCseed * seed = new AliTPCseed();
2863 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2864 Double_t cs=cos(alpha), sn=sin(alpha);
2866 // Double_t x1 =fOuterSec->GetX(i1);
2867 //Double_t xx2=fOuterSec->GetX(i2);
2869 Double_t x1 =GetXrow(i1);
2870 Double_t xx2=GetXrow(i2);
2872 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2874 Int_t imiddle = (i2+i1)/2; //middle pad row index
2875 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2876 const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2880 const AliTPCRow& kr1=GetRow(ns,i1);
2881 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2882 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2885 // change cut on curvature if it can't reach this layer
2886 // maximal curvature set to reach it
2887 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2888 if (dvertexmax*0.5*cuts[0]>0.85){
2889 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2891 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2894 if (deltay>0) ddsec = 0;
2895 // loop over clusters
2896 for (Int_t is=0; is < kr1; is++) {
2898 if (kr1[is]->IsUsed(10)) continue;
2899 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2900 //if (TMath::Abs(y1)>ymax) continue;
2902 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2904 // find possible directions
2905 Float_t anglez = (z1-z3)/(x1-x3);
2906 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2909 //find rotation angles relative to line given by vertex and point 1
2910 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2911 Double_t dvertex = TMath::Sqrt(dvertex2);
2912 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2913 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2916 // loop over 2 sectors
2922 Double_t dddz1=0; // direction of delta inclination in z axis
2929 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2930 Int_t sec2 = sec + dsec;
2932 // AliTPCRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2933 //AliTPCRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2934 AliTPCRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2935 AliTPCRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2936 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2937 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2939 // rotation angles to p1-p3
2940 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2941 Double_t x2, y2, z2;
2943 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2946 Double_t dxx0 = (xx2-x3)*cs13r;
2947 Double_t dyy0 = (xx2-x3)*sn13r;
2948 for (Int_t js=index1; js < index2; js++) {
2949 const AliTPCclusterMI *kcl = kr2[js];
2950 if (kcl->IsUsed(10)) continue;
2952 //calcutate parameters
2954 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2956 if (TMath::Abs(yy0)<0.000001) continue;
2957 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2958 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2959 Double_t r02 = (0.25+y0*y0)*dvertex2;
2960 //curvature (radius) cut
2961 if (r02<r2min) continue;
2965 Double_t c0 = 1/TMath::Sqrt(r02);
2969 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2970 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2971 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2972 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2975 Double_t z0 = kcl->GetZ();
2976 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2977 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2980 Double_t dip = (z1-z0)*c0/dfi1;
2981 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2992 x2= xx2*cs-y2*sn*dsec;
2993 y2=+xx2*sn*dsec+y2*cs;
3003 // do we have cluster at the middle ?
3005 GetProlongation(x1,xm,x,ym,zm);
3007 AliTPCclusterMI * cm=0;
3008 if (TMath::Abs(ym)-ymaxm<0){
3009 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3010 if ((!cm) || (cm->IsUsed(10))) {
3015 // rotate y1 to system 0
3016 // get state vector in rotated system
3017 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3018 Double_t xr2 = x0*cs+yr1*sn*dsec;
3019 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3021 GetProlongation(xx2,xm,xr,ym,zm);
3022 if (TMath::Abs(ym)-ymaxm<0){
3023 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3024 if ((!cm) || (cm->IsUsed(10))) {
3034 dym = ym - cm->GetY();
3035 dzm = zm - cm->GetZ();
3042 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3043 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3044 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3045 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3046 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3048 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3049 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3050 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3051 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3052 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3053 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3055 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3056 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3057 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3058 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3062 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3063 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3064 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3065 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3066 c[13]=f30*sy1*f40+f32*sy2*f42;
3067 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3069 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3071 UInt_t index=kr1.GetIndex(is);
3072 seed->~AliTPCseed(); // this does not set the pointer to 0...
3073 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3075 track->SetIsSeeding(kTRUE);
3076 track->SetSeed1(i1);
3077 track->SetSeed2(i2);
3078 track->SetSeedType(3);
3082 FollowProlongation(*track, (i1+i2)/2,1);
3083 Int_t foundable,found,shared;
3084 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3085 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3087 seed->~AliTPCseed();
3093 FollowProlongation(*track, i2,1);
3097 track->SetBConstrain(1);
3098 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3099 track->SetLastPoint(i1); // first cluster in track position
3100 track->SetFirstPoint(track->GetLastPoint());
3102 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3103 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3104 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3106 seed->~AliTPCseed();
3110 // Z VERTEX CONDITION
3111 Double_t zv, bz=GetBz();
3112 if ( !track->GetZAt(0.,bz,zv) ) continue;
3113 if (TMath::Abs(zv-z3)>cuts[2]) {
3114 FollowProlongation(*track, TMath::Max(i2-20,0));
3115 if ( !track->GetZAt(0.,bz,zv) ) continue;
3116 if (TMath::Abs(zv-z3)>cuts[2]){
3117 FollowProlongation(*track, TMath::Max(i2-40,0));
3118 if ( !track->GetZAt(0.,bz,zv) ) continue;
3119 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3120 // make seed without constrain
3121 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3122 FollowProlongation(*track2, i2,1);
3123 track2->SetBConstrain(kFALSE);
3124 track2->SetSeedType(1);
3125 arr->AddLast(track2);
3127 seed->~AliTPCseed();
3132 seed->~AliTPCseed();
3139 track->SetSeedType(0);
3140 arr->AddLast(track);
3141 seed = new AliTPCseed;
3143 // don't consider other combinations
3144 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3150 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3156 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3161 //-----------------------------------------------------------------
3162 // This function creates track seeds.
3163 //-----------------------------------------------------------------
3164 // cuts[0] - fP4 cut
3165 // cuts[1] - tan(phi) cut
3166 // cuts[2] - zvertex cut
3167 // cuts[3] - fP3 cut
3177 Double_t x[5], c[15];
3179 // make temporary seed
3180 AliTPCseed * seed = new AliTPCseed;
3181 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3182 // Double_t cs=cos(alpha), sn=sin(alpha);
3187 Double_t x1 = GetXrow(i1-1);
3188 const AliTPCRow& kr1=GetRow(sec,i1-1);
3189 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3191 Double_t x1p = GetXrow(i1);
3192 const AliTPCRow& kr1p=GetRow(sec,i1);
3194 Double_t x1m = GetXrow(i1-2);
3195 const AliTPCRow& kr1m=GetRow(sec,i1-2);
3198 //last 3 padrow for seeding
3199 AliTPCRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3200 Double_t x3 = GetXrow(i1-7);
3201 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3203 AliTPCRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3204 Double_t x3p = GetXrow(i1-6);
3206 AliTPCRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3207 Double_t x3m = GetXrow(i1-8);
3212 Int_t im = i1-4; //middle pad row index
3213 Double_t xm = GetXrow(im); // radius of middle pad-row
3214 const AliTPCRow& krm=GetRow(sec,im); //middle pad -row
3215 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3218 Double_t deltax = x1-x3;
3219 Double_t dymax = deltax*cuts[1];
3220 Double_t dzmax = deltax*cuts[3];
3222 // loop over clusters
3223 for (Int_t is=0; is < kr1; is++) {
3225 if (kr1[is]->IsUsed(10)) continue;
3226 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3228 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3230 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3231 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3237 for (Int_t js=index1; js < index2; js++) {
3238 const AliTPCclusterMI *kcl = kr3[js];
3239 if (kcl->IsUsed(10)) continue;
3241 // apply angular cuts
3242 if (TMath::Abs(y1-y3)>dymax) continue;
3245 if (TMath::Abs(z1-z3)>dzmax) continue;
3247 Double_t angley = (y1-y3)/(x1-x3);
3248 Double_t anglez = (z1-z3)/(x1-x3);
3250 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3251 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3253 Double_t yyym = angley*(xm-x1)+y1;
3254 Double_t zzzm = anglez*(xm-x1)+z1;
3256 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3258 if (kcm->IsUsed(10)) continue;
3260 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3261 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3268 // look around first
3269 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3275 if (kc1m->IsUsed(10)) used++;
3277 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3283 if (kc1p->IsUsed(10)) used++;
3285 if (used>1) continue;
3286 if (found<1) continue;
3290 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3296 if (kc3m->IsUsed(10)) used++;
3300 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3306 if (kc3p->IsUsed(10)) used++;
3310 if (used>1) continue;
3311 if (found<3) continue;
3321 x[4]=F1(x1,y1,x2,y2,x3,y3);
3322 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3325 x[2]=F2(x1,y1,x2,y2,x3,y3);
3328 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3329 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3333 Double_t sy1=0.1, sz1=0.1;
3334 Double_t sy2=0.1, sz2=0.1;
3335 Double_t sy3=0.1, sy=0.1, sz=0.1;
3337 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3338 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3339 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3340 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3341 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3342 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3344 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3345 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3346 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3347 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3351 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3352 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3353 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3354 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3355 c[13]=f30*sy1*f40+f32*sy2*f42;
3356 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3358 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3360 UInt_t index=kr1.GetIndex(is);
3361 seed->~AliTPCseed();
3362 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3364 track->SetIsSeeding(kTRUE);
3367 FollowProlongation(*track, i1-7,1);
3368 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3369 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3371 seed->~AliTPCseed();
3377 FollowProlongation(*track, i2,1);
3378 track->SetBConstrain(0);
3379 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3380 track->SetFirstPoint(track->GetLastPoint());
3382 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3383 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3384 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3386 seed->~AliTPCseed();
3391 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3392 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3393 FollowProlongation(*track2, i2,1);
3394 track2->SetBConstrain(kFALSE);
3395 track2->SetSeedType(4);
3396 arr->AddLast(track2);
3398 seed->~AliTPCseed();
3402 //arr->AddLast(track);
3403 //seed = new AliTPCseed;
3409 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3415 //_____________________________________________________________________________
3416 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3417 Float_t deltay, Bool_t /*bconstrain*/) {
3418 //-----------------------------------------------------------------
3419 // This function creates track seeds - without vertex constraint
3420 //-----------------------------------------------------------------
3421 // cuts[0] - fP4 cut - not applied
3422 // cuts[1] - tan(phi) cut
3423 // cuts[2] - zvertex cut - not applied
3424 // cuts[3] - fP3 cut
3434 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3435 // Double_t cs=cos(alpha), sn=sin(alpha);
3436 Int_t row0 = (i1+i2)/2;
3437 Int_t drow = (i1-i2)/2;
3438 const AliTPCRow& kr0=fSectors[sec][row0];
3441 AliTPCpolyTrack polytrack;
3442 Int_t nclusters=fSectors[sec][row0];
3443 AliTPCseed * seed = new AliTPCseed;
3448 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3450 Int_t nfoundable =0;
3451 for (Int_t iter =1; iter<2; iter++){ //iterations
3452 const AliTPCRow& krm=fSectors[sec][row0-iter];
3453 const AliTPCRow& krp=fSectors[sec][row0+iter];
3454 const AliTPCclusterMI * cl= kr0[is];
3456 if (cl->IsUsed(10)) {
3462 Double_t x = kr0.GetX();
3463 // Initialization of the polytrack
3468 Double_t y0= cl->GetY();
3469 Double_t z0= cl->GetZ();
3473 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3474 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3476 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3477 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3478 polytrack.AddPoint(x,y0,z0,erry, errz);
3481 if (cl->IsUsed(10)) sumused++;
3484 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3485 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3488 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3489 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3490 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3491 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3492 if (cl1->IsUsed(10)) sumused++;
3493 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3497 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3499 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3500 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3501 if (cl2->IsUsed(10)) sumused++;
3502 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3505 if (sumused>0) continue;
3507 polytrack.UpdateParameters();
3513 nfoundable = polytrack.GetN();
3514 nfound = nfoundable;
3516 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3517 Float_t maxdist = 0.8*(1.+3./(ddrow));
3518 for (Int_t delta = -1;delta<=1;delta+=2){
3519 Int_t row = row0+ddrow*delta;
3520 kr = &(fSectors[sec][row]);
3521 Double_t xn = kr->GetX();
3522 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3523 polytrack.GetFitPoint(xn,yn,zn);
3524 if (TMath::Abs(yn)>ymax) continue;
3526 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3528 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3531 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3532 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3533 if (cln->IsUsed(10)) {
3534 // printf("used\n");
3542 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3547 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3548 polytrack.UpdateParameters();
3551 if ( (sumused>3) || (sumused>0.5*nfound)) {
3552 //printf("sumused %d\n",sumused);
3557 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3558 AliTPCpolyTrack track2;
3560 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3561 if (track2.GetN()<0.5*nfoundable) continue;
3564 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3566 // test seed with and without constrain
3567 for (Int_t constrain=0; constrain<=0;constrain++){
3568 // add polytrack candidate
3570 Double_t x[5], c[15];
3571 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3572 track2.GetBoundaries(x3,x1);
3574 track2.GetFitPoint(x1,y1,z1);
3575 track2.GetFitPoint(x2,y2,z2);
3576 track2.GetFitPoint(x3,y3,z3);
3578 //is track pointing to the vertex ?
3581 polytrack.GetFitPoint(x0,y0,z0);
3594 x[4]=F1(x1,y1,x2,y2,x3,y3);
3596 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3597 x[2]=F2(x1,y1,x2,y2,x3,y3);
3599 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3600 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3601 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3602 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3605 Double_t sy =0.1, sz =0.1;
3606 Double_t sy1=0.02, sz1=0.02;
3607 Double_t sy2=0.02, sz2=0.02;
3611 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3614 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3615 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3616 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3617 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3618 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3619 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3621 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3622 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3623 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3624 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3629 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3630 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3631 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3632 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3633 c[13]=f30*sy1*f40+f32*sy2*f42;
3634 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3636 //Int_t row1 = fSectors->GetRowNumber(x1);
3637 Int_t row1 = GetRowNumber(x1);
3641 seed->~AliTPCseed();
3642 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3643 track->SetIsSeeding(kTRUE);
3644 Int_t rc=FollowProlongation(*track, i2);
3645 if (constrain) track->SetBConstrain(1);
3647 track->SetBConstrain(0);
3648 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3649 track->SetFirstPoint(track->GetLastPoint());
3651 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3652 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3653 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3656 seed->~AliTPCseed();
3659 arr->AddLast(track);
3660 seed = new AliTPCseed;
3664 } // if accepted seed
3667 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3673 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3677 //reseed using track points
3678 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3679 Int_t p1 = int(r1*track->GetNumberOfClusters());
3680 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3682 Double_t x0[3],x1[3],x2[3];
3683 for (Int_t i=0;i<3;i++){
3689 // find track position at given ratio of the length
3690 Int_t sec0=0, sec1=0, sec2=0;
3693 for (Int_t i=0;i<160;i++){
3694 if (track->GetClusterPointer(i)){
3696 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3697 if ( (index<p0) || x0[0]<0 ){
3698 if (trpoint->GetX()>1){
3699 clindex = track->GetClusterIndex2(i);
3701 x0[0] = trpoint->GetX();
3702 x0[1] = trpoint->GetY();
3703 x0[2] = trpoint->GetZ();
3704 sec0 = ((clindex&0xff000000)>>24)%18;
3709 if ( (index<p1) &&(trpoint->GetX()>1)){
3710 clindex = track->GetClusterIndex2(i);
3712 x1[0] = trpoint->GetX();
3713 x1[1] = trpoint->GetY();
3714 x1[2] = trpoint->GetZ();
3715 sec1 = ((clindex&0xff000000)>>24)%18;
3718 if ( (index<p2) &&(trpoint->GetX()>1)){
3719 clindex = track->GetClusterIndex2(i);
3721 x2[0] = trpoint->GetX();
3722 x2[1] = trpoint->GetY();
3723 x2[2] = trpoint->GetZ();
3724 sec2 = ((clindex&0xff000000)>>24)%18;
3731 Double_t alpha, cs,sn, xx2,yy2;
3733 alpha = (sec1-sec2)*fSectors->GetAlpha();
3734 cs = TMath::Cos(alpha);
3735 sn = TMath::Sin(alpha);
3736 xx2= x1[0]*cs-x1[1]*sn;
3737 yy2= x1[0]*sn+x1[1]*cs;
3741 alpha = (sec0-sec2)*fSectors->GetAlpha();
3742 cs = TMath::Cos(alpha);
3743 sn = TMath::Sin(alpha);
3744 xx2= x0[0]*cs-x0[1]*sn;
3745 yy2= x0[0]*sn+x0[1]*cs;
3751 Double_t x[5],c[15];
3755 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3756 // if (x[4]>1) return 0;
3757 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3758 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3759 //if (TMath::Abs(x[3]) > 2.2) return 0;
3760 //if (TMath::Abs(x[2]) > 1.99) return 0;
3762 Double_t sy =0.1, sz =0.1;
3764 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3765 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3766 Double_t sy3=0.01+track->GetSigmaY2();
3768 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3769 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3770 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3771 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3772 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3773 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3775 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3776 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3777 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3778 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3783 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3784 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3785 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3786 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3787 c[13]=f30*sy1*f40+f32*sy2*f42;
3788 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3790 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3791 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3792 // Double_t y0,z0,y1,z1, y2,z2;
3793 //seed->GetProlongation(x0[0],y0,z0);
3794 // seed->GetProlongation(x1[0],y1,z1);
3795 //seed->GetProlongation(x2[0],y2,z2);
3797 seed->SetLastPoint(pp2);
3798 seed->SetFirstPoint(pp2);
3805 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3809 //reseed using founded clusters
3811 // Find the number of clusters
3812 Int_t nclusters = 0;
3813 for (Int_t irow=0;irow<160;irow++){
3814 if (track->GetClusterIndex(irow)>0) nclusters++;
3818 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3819 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3820 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3824 Int_t row[3],sec[3]={0,0,0};
3826 // find track row position at given ratio of the length
3828 for (Int_t irow=0;irow<160;irow++){
3829 if (track->GetClusterIndex2(irow)<0) continue;
3831 for (Int_t ipoint=0;ipoint<3;ipoint++){
3832 if (index<=ipos[ipoint]) row[ipoint] = irow;
3836 //Get cluster and sector position
3837 for (Int_t ipoint=0;ipoint<3;ipoint++){
3838 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3839 AliTPCclusterMI * cl = GetClusterMI(clindex);
3842 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3845 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3846 xyz[ipoint][0] = GetXrow(row[ipoint]);
3847 xyz[ipoint][1] = cl->GetY();
3848 xyz[ipoint][2] = cl->GetZ();
3852 // Calculate seed state vector and covariance matrix
3854 Double_t alpha, cs,sn, xx2,yy2;
3856 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3857 cs = TMath::Cos(alpha);
3858 sn = TMath::Sin(alpha);
3859 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3860 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3864 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3865 cs = TMath::Cos(alpha);
3866 sn = TMath::Sin(alpha);
3867 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3868 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3874 Double_t x[5],c[15];
3878 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3879 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3880 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3882 Double_t sy =0.1, sz =0.1;
3884 Double_t sy1=0.2, sz1=0.2;
3885 Double_t sy2=0.2, sz2=0.2;
3888 Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
3889 Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
3890 Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
3891 Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
3892 Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
3893 Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
3895 Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
3896 Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
3897 Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
3898 Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
3903 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3904 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3905 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3906 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3907 c[13]=f30*sy1*f40+f32*sy2*f42;
3908 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3910 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3911 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3912 seed->SetLastPoint(row[2]);
3913 seed->SetFirstPoint(row[2]);
3918 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3922 //reseed using founded clusters
3925 Int_t row[3]={0,0,0};
3926 Int_t sec[3]={0,0,0};
3928 // forward direction
3930 for (Int_t irow=r0;irow<160;irow++){
3931 if (track->GetClusterIndex(irow)>0){
3936 for (Int_t irow=160;irow>r0;irow--){
3937 if (track->GetClusterIndex(irow)>0){
3942 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3943 if (track->GetClusterIndex(irow)>0){
3951 for (Int_t irow=0;irow<r0;irow++){
3952 if (track->GetClusterIndex(irow)>0){
3957 for (Int_t irow=r0;irow>0;irow--){
3958 if (track->GetClusterIndex(irow)>0){
3963 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3964 if (track->GetClusterIndex(irow)>0){
3971 if ((row[2]-row[0])<20) return 0;
3972 if (row[1]==0) return 0;
3975 //Get cluster and sector position
3976 for (Int_t ipoint=0;ipoint<3;ipoint++){
3977 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3978 AliTPCclusterMI * cl = GetClusterMI(clindex);
3981 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3984 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3985 xyz[ipoint][0] = GetXrow(row[ipoint]);
3986 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3987 if (point&&ipoint<2){
3989 xyz[ipoint][1] = point->GetY();
3990 xyz[ipoint][2] = point->GetZ();
3993 xyz[ipoint][1] = cl->GetY();
3994 xyz[ipoint][2] = cl->GetZ();
4001 // Calculate seed state vector and covariance matrix
4003 Double_t alpha, cs,sn, xx2,yy2;
4005 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4006 cs = TMath::Cos(alpha);
4007 sn = TMath::Sin(alpha);
4008 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4009 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4013 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4014 cs = TMath::Cos(alpha);
4015 sn = TMath::Sin(alpha);
4016 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4017 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4023 Double_t x[5],c[15];
4027 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4028 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4029 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4031 Double_t sy =0.1, sz =0.1;
4033 Double_t sy1=0.2, sz1=0.2;
4034 Double_t sy2=0.2, sz2=0.2;
4037 Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
4038 Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
4039 Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
4040 Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
4041 Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
4042 Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
4044 Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
4045 Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
4046 Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
4047 Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
4052 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4053 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4054 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4055 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4056 c[13]=f30*sy1*f40+f32*sy2*f42;
4057 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4059 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4060 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4061 seed->SetLastPoint(row[2]);
4062 seed->SetFirstPoint(row[2]);
4063 for (Int_t i=row[0];i<row[2];i++){
4064 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4072 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4075 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4077 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4079 // Two reasons to have multiple find tracks
4080 // 1. Curling tracks can be find more than once
4081 // 2. Splitted tracks
4082 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4083 // b.) Edge effect on the sector boundaries
4086 // Algorithm done in 2 phases - because of CPU consumption
4087 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4089 // Algorihm for curling tracks sign:
4090 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4091 // a.) opposite sign
4092 // b.) one of the tracks - not pointing to the primary vertex -
4093 // c.) delta tan(theta)
4095 // 2 phase - calculates DCA between tracks - time consument
4100 // General cuts - for splitted tracks and for curling tracks
4102 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4104 // Curling tracks cuts
4109 Int_t nentries = array->GetEntriesFast();
4110 AliHelix *helixes = new AliHelix[nentries];
4111 Float_t *xm = new Float_t[nentries];
4112 Float_t *dz0 = new Float_t[nentries];
4113 Float_t *dz1 = new Float_t[nentries];
4119 // Find track COG in x direction - point with best defined parameters
4121 for (Int_t i=0;i<nentries;i++){
4122 AliTPCseed* track = (AliTPCseed*)array->At(i);
4123 if (!track) continue;
4124 track->SetCircular(0);
4125 new (&helixes[i]) AliHelix(*track);
4129 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4132 for (Int_t icl=0; icl<160; icl++){
4133 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4139 if (ncl>0) xm[i]/=Float_t(ncl);
4141 TTreeSRedirector &cstream = *fDebugStreamer;
4143 for (Int_t i0=0;i0<nentries;i0++){
4144 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4145 if (!track0) continue;
4146 Float_t xc0 = helixes[i0].GetHelix(6);
4147 Float_t yc0 = helixes[i0].GetHelix(7);
4148 Float_t r0 = helixes[i0].GetHelix(8);
4149 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4150 Float_t fi0 = TMath::ATan2(yc0,xc0);
4152 for (Int_t i1=i0+1;i1<nentries;i1++){
4153 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4154 if (!track1) continue;
4155 Int_t lab0=track0->GetLabel();
4156 Int_t lab1=track1->GetLabel();
4157 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4159 Float_t xc1 = helixes[i1].GetHelix(6);
4160 Float_t yc1 = helixes[i1].GetHelix(7);
4161 Float_t r1 = helixes[i1].GetHelix(8);
4162 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4163 Float_t fi1 = TMath::ATan2(yc1,xc1);
4165 Float_t dfi = fi0-fi1;
4168 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4169 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4170 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4172 // if short tracks with undefined sign
4173 fi1 = -TMath::ATan2(yc1,-xc1);
4176 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4179 // debug stream to tune "fast cuts"
4181 Double_t dist[3]; // distance at X
4182 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4183 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4184 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4185 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4186 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4187 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4188 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4189 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4193 for (Int_t icl=0; icl<160; icl++){
4194 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4195 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4198 if (cl0==cl1) sums++;
4206 "Tr0.="<<track0<< // seed0
4207 "Tr1.="<<track1<< // seed1
4208 "h0.="<<&helixes[i0]<<
4209 "h1.="<<&helixes[i1]<<
4211 "sum="<<sum<< //the sum of rows with cl in both
4212 "sums="<<sums<< //the sum of shared clusters
4213 "xm0="<<xm[i0]<< // the center of track
4214 "xm1="<<xm[i1]<< // the x center of track
4215 // General cut variables
4216 "dfi="<<dfi<< // distance in fi angle
4217 "dtheta="<<dtheta<< // distance int theta angle
4223 "dist0="<<dist[0]<< //distance x
4224 "dist1="<<dist[1]<< //distance y
4225 "dist2="<<dist[2]<< //distance z
4226 "mdist0="<<mdist[0]<< //distance x
4227 "mdist1="<<mdist[1]<< //distance y
4228 "mdist2="<<mdist[2]<< //distance z
4241 if (AliTPCReconstructor::StreamLevel()>1) {
4242 AliInfo("Time for curling tracks removal DEBUGGING MC");
4248 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4252 // Two reasons to have multiple find tracks
4253 // 1. Curling tracks can be find more than once
4254 // 2. Splitted tracks
4255 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4256 // b.) Edge effect on the sector boundaries
4258 // This function tries to find tracks closed in the parametric space
4260 // cut logic if distance is bigger than cut continue - Do Nothing
4261 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4262 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4263 const Float_t kdelta = 40.; //delta r to calculate track distance
4265 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4266 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4269 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4270 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4275 Int_t nentries = array->GetEntriesFast();
4276 AliHelix *helixes = new AliHelix[nentries];
4277 Float_t *xm = new Float_t[nentries];
4283 //Sort tracks according quality
4285 Int_t nseed = array->GetEntriesFast();
4286 Float_t * quality = new Float_t[nseed];
4287 Int_t * indexes = new Int_t[nseed];
4288 for (Int_t i=0; i<nseed; i++) {
4289 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4294 pt->UpdatePoints(); //select first last max dens points
4295 Float_t * points = pt->GetPoints();
4296 if (points[3]<0.8) quality[i] =-1;
4297 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4298 //prefer high momenta tracks if overlaps
4299 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4301 TMath::Sort(nseed,quality,indexes);
4305 // Find track COG in x direction - point with best defined parameters
4307 for (Int_t i=0;i<nentries;i++){
4308 AliTPCseed* track = (AliTPCseed*)array->At(i);
4309 if (!track) continue;
4310 track->SetCircular(0);
4311 new (&helixes[i]) AliHelix(*track);
4314 for (Int_t icl=0; icl<160; icl++){
4315 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4321 if (ncl>0) xm[i]/=Float_t(ncl);
4323 TTreeSRedirector &cstream = *fDebugStreamer;
4325 for (Int_t is0=0;is0<nentries;is0++){
4326 Int_t i0 = indexes[is0];
4327 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4328 if (!track0) continue;
4329 if (track0->GetKinkIndexes()[0]!=0) continue;
4330 Float_t xc0 = helixes[i0].GetHelix(6);
4331 Float_t yc0 = helixes[i0].GetHelix(7);
4332 Float_t fi0 = TMath::ATan2(yc0,xc0);
4334 for (Int_t is1=is0+1;is1<nentries;is1++){
4335 Int_t i1 = indexes[is1];
4336 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4337 if (!track1) continue;
4339 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4340 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4341 if (track1->GetKinkIndexes()[0]!=0) continue;
4343 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4344 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4346 Float_t xc1 = helixes[i1].GetHelix(6);
4347 Float_t yc1 = helixes[i1].GetHelix(7);
4348 Float_t fi1 = TMath::ATan2(yc1,xc1);
4350 Float_t dfi = fi0-fi1;
4351 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4352 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4353 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4355 // if short tracks with undefined sign
4356 fi1 = -TMath::ATan2(yc1,-xc1);
4359 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4366 for (Int_t icl=0; icl<160; icl++){
4367 Int_t index0=track0->GetClusterIndex2(icl);
4368 Int_t index1=track1->GetClusterIndex2(icl);
4369 Bool_t used0 = (index0>0 && !(index0&0x8000));
4370 Bool_t used1 = (index1>0 && !(index1&0x8000));
4372 if (used0) sum0++; // used cluster0
4373 if (used1) sum1++; // used clusters1
4374 if (used0&&used1) sum++;
4375 if (index0==index1 && used0 && used1) sums++;
4379 if (sums<10) continue;
4380 if (sum<40) continue;
4381 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4383 Double_t dist[5][4]; // distance at X
4384 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4388 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4389 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4390 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4391 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4393 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4394 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4395 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4396 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4398 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4399 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4400 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4403 Int_t lab0=track0->GetLabel();
4404 Int_t lab1=track1->GetLabel();
4405 cstream<<"Splitted"<<
4409 "Tr0.="<<track0<< // seed0
4410 "Tr1.="<<track1<< // seed1
4411 "h0.="<<&helixes[i0]<<
4412 "h1.="<<&helixes[i1]<<
4414 "sum="<<sum<< //the sum of rows with cl in both
4415 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4416 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4417 "sums="<<sums<< //the sum of shared clusters
4418 "xm0="<<xm[i0]<< // the center of track
4419 "xm1="<<xm[i1]<< // the x center of track
4420 // General cut variables
4421 "dfi="<<dfi<< // distance in fi angle
4422 "dtheta="<<dtheta<< // distance int theta angle
4425 "dist0="<<dist[4][0]<< //distance x
4426 "dist1="<<dist[4][1]<< //distance y
4427 "dist2="<<dist[4][2]<< //distance z
4428 "mdist0="<<mdist[0]<< //distance x
4429 "mdist1="<<mdist[1]<< //distance y
4430 "mdist2="<<mdist[2]<< //distance z
4433 delete array->RemoveAt(i1);
4438 AliInfo("Time for splitted tracks removal");
4444 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4447 // find Curling tracks
4448 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4451 // Algorithm done in 2 phases - because of CPU consumption
4452 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4453 // see detal in MC part what can be used to cut
4457 const Float_t kMaxC = 400; // maximal curvature to of the track
4458 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4459 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4460 const Float_t kPtRatio = 0.3; // ratio between pt
4461 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4464 // Curling tracks cuts
4467 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4468 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4469 const Float_t kMinAngle = 2.9; // angle between tracks
4470 const Float_t kMaxDist = 5; // biggest distance
4472 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4475 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4476 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4477 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4478 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4479 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4481 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4482 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4484 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4485 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4487 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4493 Int_t nentries = array->GetEntriesFast();
4494 AliHelix *helixes = new AliHelix[nentries];
4495 for (Int_t i=0;i<nentries;i++){
4496 AliTPCseed* track = (AliTPCseed*)array->At(i);
4497 if (!track) continue;
4498 track->SetCircular(0);
4499 new (&helixes[i]) AliHelix(*track);
4505 Double_t phase[2][2],radius[2];
4509 TTreeSRedirector &cstream = *fDebugStreamer;
4511 for (Int_t i0=0;i0<nentries;i0++){
4512 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4513 if (!track0) continue;
4514 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4515 Float_t xc0 = helixes[i0].GetHelix(6);
4516 Float_t yc0 = helixes[i0].GetHelix(7);
4517 Float_t r0 = helixes[i0].GetHelix(8);
4518 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4519 Float_t fi0 = TMath::ATan2(yc0,xc0);
4521 for (Int_t i1=i0+1;i1<nentries;i1++){
4522 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4523 if (!track1) continue;
4524 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4525 Float_t xc1 = helixes[i1].GetHelix(6);
4526 Float_t yc1 = helixes[i1].GetHelix(7);
4527 Float_t r1 = helixes[i1].GetHelix(8);
4528 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4529 Float_t fi1 = TMath::ATan2(yc1,xc1);
4531 Float_t dfi = fi0-fi1;
4534 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4535 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4536 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4540 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4541 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4542 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4543 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4544 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4546 Float_t pt0 = track0->GetSignedPt();
4547 Float_t pt1 = track1->GetSignedPt();
4548 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4549 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4550 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4551 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4554 // Now find closest approach
4558 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4559 if (npoints==0) continue;
4560 helixes[i0].GetClosestPhases(helixes[i1], phase);
4564 Double_t hangles[3];
4565 helixes[i0].Evaluate(phase[0][0],xyz0);
4566 helixes[i1].Evaluate(phase[0][1],xyz1);
4568 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4569 Double_t deltah[2],deltabest;
4570 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4574 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4576 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4577 if (deltah[1]<deltah[0]) ibest=1;
4579 deltabest = TMath::Sqrt(deltah[ibest]);
4580 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4581 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4582 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4583 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4585 if (deltabest>kMaxDist) continue;
4586 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4587 Bool_t sign =kFALSE;
4588 if (hangles[2]>kMinAngle) sign =kTRUE;
4591 // circular[i0] = kTRUE;
4592 // circular[i1] = kTRUE;
4593 if (track0->OneOverPt()<track1->OneOverPt()){
4594 track0->SetCircular(track0->GetCircular()+1);
4595 track1->SetCircular(track1->GetCircular()+2);
4598 track1->SetCircular(track1->GetCircular()+1);
4599 track0->SetCircular(track0->GetCircular()+2);
4602 if (AliTPCReconstructor::StreamLevel()>1){
4604 //debug stream to tune "fine" cuts
4605 Int_t lab0=track0->GetLabel();
4606 Int_t lab1=track1->GetLabel();
4607 cstream<<"Curling2"<<
4623 "npoints="<<npoints<<
4624 "hangles0="<<hangles[0]<<
4625 "hangles1="<<hangles[1]<<
4626 "hangles2="<<hangles[2]<<
4629 "radius="<<radiusbest<<
4630 "deltabest="<<deltabest<<
4631 "phase0="<<phase[ibest][0]<<
4632 "phase1="<<phase[ibest][1]<<
4640 if (AliTPCReconstructor::StreamLevel()>1) {
4641 AliInfo("Time for curling tracks removal");
4650 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4657 TObjArray *kinks= new TObjArray(10000);
4658 // TObjArray *v0s= new TObjArray(10000);
4659 Int_t nentries = array->GetEntriesFast();
4660 AliHelix *helixes = new AliHelix[nentries];
4661 Int_t *sign = new Int_t[nentries];
4662 Int_t *nclusters = new Int_t[nentries];
4663 Float_t *alpha = new Float_t[nentries];
4664 AliKink *kink = new AliKink();
4665 Int_t * usage = new Int_t[nentries];
4666 Float_t *zm = new Float_t[nentries];
4667 Float_t *z0 = new Float_t[nentries];
4668 Float_t *fim = new Float_t[nentries];
4669 Float_t *shared = new Float_t[nentries];
4670 Bool_t *circular = new Bool_t[nentries];
4671 Float_t *dca = new Float_t[nentries];
4672 //const AliESDVertex * primvertex = esd->GetVertex();
4674 // nentries = array->GetEntriesFast();
4679 for (Int_t i=0;i<nentries;i++){
4682 AliTPCseed* track = (AliTPCseed*)array->At(i);
4683 if (!track) continue;
4684 track->SetCircular(0);
4686 track->UpdatePoints();
4687 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4689 nclusters[i]=track->GetNumberOfClusters();
4690 alpha[i] = track->GetAlpha();
4691 new (&helixes[i]) AliHelix(*track);
4693 helixes[i].Evaluate(0,xyz);
4694 sign[i] = (track->GetC()>0) ? -1:1;
4697 if (track->GetProlongation(x,y,z)){
4699 fim[i] = alpha[i]+TMath::ATan2(y,x);
4702 zm[i] = track->GetZ();
4706 circular[i]= kFALSE;
4707 if (track->GetProlongation(0,y,z)) z0[i] = z;
4708 dca[i] = track->GetD(0,0);
4714 Int_t ncandidates =0;
4717 Double_t phase[2][2],radius[2];
4720 // Find circling track
4721 TTreeSRedirector &cstream = *fDebugStreamer;
4723 for (Int_t i0=0;i0<nentries;i0++){
4724 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4725 if (!track0) continue;
4726 if (track0->GetNumberOfClusters()<40) continue;
4727 if (TMath::Abs(1./track0->GetC())>200) continue;
4728 for (Int_t i1=i0+1;i1<nentries;i1++){
4729 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4730 if (!track1) continue;
4731 if (track1->GetNumberOfClusters()<40) continue;
4732 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4733 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4734 if (TMath::Abs(1./track1->GetC())>200) continue;
4735 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4736 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4737 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4738 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4739 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4741 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4742 if (mindcar<5) continue;
4743 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4744 if (mindcaz<5) continue;
4745 if (mindcar+mindcaz<20) continue;
4748 Float_t xc0 = helixes[i0].GetHelix(6);
4749 Float_t yc0 = helixes[i0].GetHelix(7);
4750 Float_t r0 = helixes[i0].GetHelix(8);
4751 Float_t xc1 = helixes[i1].GetHelix(6);
4752 Float_t yc1 = helixes[i1].GetHelix(7);
4753 Float_t r1 = helixes[i1].GetHelix(8);
4755 Float_t rmean = (r0+r1)*0.5;
4756 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4757 //if (delta>30) continue;
4758 if (delta>rmean*0.25) continue;
4759 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4761 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4762 if (npoints==0) continue;
4763 helixes[i0].GetClosestPhases(helixes[i1], phase);
4767 Double_t hangles[3];
4768 helixes[i0].Evaluate(phase[0][0],xyz0);
4769 helixes[i1].Evaluate(phase[0][1],xyz1);
4771 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4772 Double_t deltah[2],deltabest;
4773 if (hangles[2]<2.8) continue;
4776 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4778 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4779 if (deltah[1]<deltah[0]) ibest=1;
4781 deltabest = TMath::Sqrt(deltah[ibest]);
4782 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4783 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4784 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4785 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4787 if (deltabest>6) continue;
4788 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4789 Bool_t sign =kFALSE;
4790 if (hangles[2]>3.06) sign =kTRUE;
4793 circular[i0] = kTRUE;
4794 circular[i1] = kTRUE;
4795 if (track0->OneOverPt()<track1->OneOverPt()){
4796 track0->SetCircular(track0->GetCircular()+1);
4797 track1->SetCircular(track1->GetCircular()+2);
4800 track1->SetCircular(track1->GetCircular()+1);
4801 track0->SetCircular(track0->GetCircular()+2);
4804 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4806 Int_t lab0=track0->GetLabel();
4807 Int_t lab1=track1->GetLabel();
4808 cstream<<"Curling"<<
4815 "mindcar="<<mindcar<<
4816 "mindcaz="<<mindcaz<<
4819 "npoints="<<npoints<<
4820 "hangles0="<<hangles[0]<<
4821 "hangles2="<<hangles[2]<<
4826 "radius="<<radiusbest<<
4827 "deltabest="<<deltabest<<
4828 "phase0="<<phase[ibest][0]<<
4829 "phase1="<<phase[ibest][1]<<
4839 for (Int_t i =0;i<nentries;i++){
4840 if (sign[i]==0) continue;
4841 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4848 Double_t cradius0 = 40*40;
4849 Double_t cradius1 = 270*270;
4852 Double_t cdist3=0.55;
4853 for (Int_t j =i+1;j<nentries;j++){
4855 if (sign[j]*sign[i]<1) continue;
4856 if ( (nclusters[i]+nclusters[j])>200) continue;
4857 if ( (nclusters[i]+nclusters[j])<80) continue;
4858 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4859 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4860 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4861 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4862 if (npoints<1) continue;
4865 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4868 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4871 Double_t delta1=10000,delta2=10000;
4872 // cuts on the intersection radius
4873 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4874 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4875 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4877 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4878 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4879 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4882 Double_t distance1 = TMath::Min(delta1,delta2);
4883 if (distance1>cdist1) continue; // cut on DCA linear approximation
4885 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4886 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4887 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4888 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4891 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4892 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4893 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4895 distance1 = TMath::Min(delta1,delta2);
4898 rkink = TMath::Sqrt(radius[0]);
4901 rkink = TMath::Sqrt(radius[1]);
4903 if (distance1>cdist2) continue;
4906 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4909 Int_t row0 = GetRowNumber(rkink);
4910 if (row0<10) continue;
4911 if (row0>150) continue;
4914 Float_t dens00=-1,dens01=-1;
4915 Float_t dens10=-1,dens11=-1;
4917 Int_t found,foundable,shared;
4918 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4919 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4920 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4921 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4923 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4924 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4925 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4926 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4928 if (dens00<dens10 && dens01<dens11) continue;
4929 if (dens00>dens10 && dens01>dens11) continue;
4930 if (TMath::Max(dens00,dens10)<0.1) continue;
4931 if (TMath::Max(dens01,dens11)<0.3) continue;
4933 if (TMath::Min(dens00,dens10)>0.6) continue;
4934 if (TMath::Min(dens01,dens11)>0.6) continue;
4937 AliTPCseed * ktrack0, *ktrack1;
4946 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4947 AliExternalTrackParam paramm(*ktrack0);
4948 AliExternalTrackParam paramd(*ktrack1);
4949 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4952 kink->SetMother(paramm);
4953 kink->SetDaughter(paramd);
4956 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4958 fParam->Transform0to1(x,index);
4959 fParam->Transform1to2(x,index);
4960 row0 = GetRowNumber(x[0]);
4962 if (kink->GetR()<100) continue;
4963 if (kink->GetR()>240) continue;
4964 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4965 if (kink->GetDistance()>cdist3) continue;
4966 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4967 if (dird<0) continue;
4969 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4970 if (dirm<0) continue;
4971 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4972 if (mpt<0.2) continue;
4975 //for high momenta momentum not defined well in first iteration
4976 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4977 if (qt>0.35) continue;
4980 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4981 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4983 kink->SetTPCDensity(dens00,0,0);
4984 kink->SetTPCDensity(dens01,0,1);
4985 kink->SetTPCDensity(dens10,1,0);
4986 kink->SetTPCDensity(dens11,1,1);
4987 kink->SetIndex(i,0);
4988 kink->SetIndex(j,1);
4991 kink->SetTPCDensity(dens10,0,0);
4992 kink->SetTPCDensity(dens11,0,1);
4993 kink->SetTPCDensity(dens00,1,0);
4994 kink->SetTPCDensity(dens01,1,1);
4995 kink->SetIndex(j,0);
4996 kink->SetIndex(i,1);
4999 if (mpt<1||kink->GetAngle(2)>0.1){
5000 // angle and densities not defined yet
5001 if (kink->GetTPCDensityFactor()<0.8) continue;
5002 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5003 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5004 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5005 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5007 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5008 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5009 criticalangle= 3*TMath::Sqrt(criticalangle);
5010 if (criticalangle>0.02) criticalangle=0.02;
5011 if (kink->GetAngle(2)<criticalangle) continue;
5014 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5015 Float_t shapesum =0;
5017 for ( Int_t row = row0-drow; row<row0+drow;row++){
5018 if (row<0) continue;
5019 if (row>155) continue;
5020 if (ktrack0->GetClusterPointer(row)){
5021 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5022 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5025 if (ktrack1->GetClusterPointer(row)){
5026 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5027 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5032 kink->SetShapeFactor(-1.);
5035 kink->SetShapeFactor(shapesum/sum);
5037 // esd->AddKink(kink);
5038 kinks->AddLast(kink);
5044 // sort the kinks according quality - and refit them towards vertex
5046 Int_t nkinks = kinks->GetEntriesFast();
5047 Float_t *quality = new Float_t[nkinks];
5048 Int_t *indexes = new Int_t[nkinks];
5049 AliTPCseed *mothers = new AliTPCseed[nkinks];
5050 AliTPCseed *daughters = new AliTPCseed[nkinks];
5053 for (Int_t i=0;i<nkinks;i++){
5055 AliKink *kink = (AliKink*)kinks->At(i);
5057 // refit kinks towards vertex
5059 Int_t index0 = kink->GetIndex(0);
5060 Int_t index1 = kink->GetIndex(1);
5061 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5062 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5064 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5066 // Refit Kink under if too small angle
5068 if (kink->GetAngle(2)<0.05){
5069 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5070 Int_t row0 = kink->GetTPCRow0();
5071 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
5074 Int_t last = row0-drow;
5075 if (last<40) last=40;
5076 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5077 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5080 Int_t first = row0+drow;
5081 if (first>130) first=130;
5082 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5083 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5085 if (seed0 && seed1){
5086 kink->SetStatus(1,8);
5087 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
5088 row0 = GetRowNumber(kink->GetR());
5089 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5090 mothers[i] = *seed0;
5091 daughters[i] = *seed1;
5094 delete kinks->RemoveAt(i);
5095 if (seed0) delete seed0;
5096 if (seed1) delete seed1;
5099 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
5100 delete kinks->RemoveAt(i);
5101 if (seed0) delete seed0;
5102 if (seed1) delete seed1;
5110 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5112 TMath::Sort(nkinks,quality,indexes,kFALSE);
5114 //remove double find kinks
5116 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5117 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5118 if (!kink0) continue;
5120 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5121 if (!kink0) continue;
5122 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5123 if (!kink1) continue;
5124 // if not close kink continue
5125 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5126 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5127 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5129 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5130 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5131 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5132 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5133 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5142 for (Int_t i=0;i<row0;i++){
5143 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5146 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5153 for (Int_t i=row0;i<158;i++){
5154 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5157 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5163 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5164 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5165 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5166 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5167 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5168 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5170 shared[kink0->GetIndex(0)]= kTRUE;
5171 shared[kink0->GetIndex(1)]= kTRUE;
5172 delete kinks->RemoveAt(indexes[ikink0]);
5175 shared[kink1->GetIndex(0)]= kTRUE;
5176 shared[kink1->GetIndex(1)]= kTRUE;
5177 delete kinks->RemoveAt(indexes[ikink1]);
5184 for (Int_t i=0;i<nkinks;i++){
5185 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5186 if (!kink) continue;
5187 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5188 Int_t index0 = kink->GetIndex(0);
5189 Int_t index1 = kink->GetIndex(1);
5190 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5191 kink->SetMultiple(usage[index0],0);
5192 kink->SetMultiple(usage[index1],1);
5193 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5194 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5195 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5196 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5198 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5199 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5200 if (!ktrack0 || !ktrack1) continue;
5201 Int_t index = esd->AddKink(kink);
5204 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5205 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5206 *ktrack0 = mothers[indexes[i]];
5207 *ktrack1 = daughters[indexes[i]];
5211 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5212 ktrack1->SetKinkIndex(usage[index1], (index+1));
5217 // Remove tracks corresponding to shared kink's
5219 for (Int_t i=0;i<nentries;i++){
5220 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5221 if (!track0) continue;
5222 if (track0->GetKinkIndex(0)!=0) continue;
5223 if (shared[i]) delete array->RemoveAt(i);
5228 RemoveUsed2(array,0.5,0.4,30);
5230 for (Int_t i=0;i<nentries;i++){
5231 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5232 if (!track0) continue;
5233 track0->CookdEdx(0.02,0.6);
5237 for (Int_t i=0;i<nentries;i++){
5238 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5239 if (!track0) continue;
5240 if (track0->Pt()<1.4) continue;
5241 //remove double high momenta tracks - overlapped with kink candidates
5244 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5245 if (track0->GetClusterPointer(icl)!=0){
5247 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5250 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5251 delete array->RemoveAt(i);
5255 if (track0->GetKinkIndex(0)!=0) continue;
5256 if (track0->GetNumberOfClusters()<80) continue;
5258 AliTPCseed *pmother = new AliTPCseed();
5259 AliTPCseed *pdaughter = new AliTPCseed();
5260 AliKink *pkink = new AliKink;
5262 AliTPCseed & mother = *pmother;
5263 AliTPCseed & daughter = *pdaughter;
5264 AliKink & kink = *pkink;
5265 if (CheckKinkPoint(track0,mother,daughter, kink)){
5266 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5270 continue; //too short tracks
5272 if (mother.Pt()<1.4) {
5278 Int_t row0= kink.GetTPCRow0();
5279 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5286 Int_t index = esd->AddKink(&kink);
5287 mother.SetKinkIndex(0,-(index+1));
5288 daughter.SetKinkIndex(0,index+1);
5289 if (mother.GetNumberOfClusters()>50) {
5290 delete array->RemoveAt(i);
5291 array->AddAt(new AliTPCseed(mother),i);
5294 array->AddLast(new AliTPCseed(mother));
5296 array->AddLast(new AliTPCseed(daughter));
5297 for (Int_t icl=0;icl<row0;icl++) {
5298 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5301 for (Int_t icl=row0;icl<158;icl++) {
5302 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5311 delete [] daughters;
5333 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5337 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5343 TObjArray *tpcv0s = new TObjArray(100000);
5344 Int_t nentries = array->GetEntriesFast();
5345 AliHelix *helixes = new AliHelix[nentries];
5346 Int_t *sign = new Int_t[nentries];
5347 Float_t *alpha = new Float_t[nentries];
5348 Float_t *z0 = new Float_t[nentries];
5349 Float_t *dca = new Float_t[nentries];
5350 Float_t *sdcar = new Float_t[nentries];
5351 Float_t *cdcar = new Float_t[nentries];
5352 Float_t *pulldcar = new Float_t[nentries];
5353 Float_t *pulldcaz = new Float_t[nentries];
5354 Float_t *pulldca = new Float_t[nentries];
5355 Bool_t *isPrim = new Bool_t[nentries];
5356 const AliESDVertex * primvertex = esd->GetVertex();
5357 Double_t zvertex = primvertex->GetZv();
5359 // nentries = array->GetEntriesFast();
5361 for (Int_t i=0;i<nentries;i++){
5364 AliTPCseed* track = (AliTPCseed*)array->At(i);
5365 if (!track) continue;
5366 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5367 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5368 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5370 alpha[i] = track->GetAlpha();
5371 new (&helixes[i]) AliHelix(*track);
5373 helixes[i].Evaluate(0,xyz);
5374 sign[i] = (track->GetC()>0) ? -1:1;
5378 if (track->GetProlongation(0,y,z)) z0[i] = z;
5379 dca[i] = track->GetD(0,0);
5381 // dca error parrameterezation + pulls
5383 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5384 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5385 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5386 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5387 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5388 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5389 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5390 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5392 if (track->TPCrPID(4)>0.5) {
5393 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5395 if (track->TPCrPID(0)>0.4) {
5396 isPrim[i]=kFALSE; //electron no sigma cut
5403 Int_t ncandidates =0;
5406 Double_t phase[2][2],radius[2];
5412 TTreeSRedirector &cstream = *fDebugStreamer;
5413 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5415 Double_t cradius0 = 10*10;
5416 Double_t cradius1 = 200*200;
5419 Double_t cpointAngle = 0.95;
5421 Double_t delta[2]={10000,10000};
5422 for (Int_t i =0;i<nentries;i++){
5423 if (sign[i]==0) continue;
5424 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5425 if (!track0) continue;
5426 if (AliTPCReconstructor::StreamLevel()>1){
5431 "zvertex="<<zvertex<<
5432 "sdcar0="<<sdcar[i]<<
5433 "cdcar0="<<cdcar[i]<<
5434 "pulldcar0="<<pulldcar[i]<<
5435 "pulldcaz0="<<pulldcaz[i]<<
5436 "pulldca0="<<pulldca[i]<<
5437 "isPrim="<<isPrim[i]<<
5441 if (track0->GetSigned1Pt()<0) continue;
5442 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5444 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5449 for (Int_t j =0;j<nentries;j++){
5450 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5451 if (!track1) continue;
5452 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5453 if (sign[j]*sign[i]>0) continue;
5454 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5455 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5458 // DCA to prim vertex cut
5464 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5465 if (npoints<1) continue;
5469 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5470 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5471 if (delta[0]>cdist1) continue;
5474 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5475 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5476 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5477 if (delta[1]<delta[0]) iclosest=1;
5478 if (delta[iclosest]>cdist1) continue;
5480 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5481 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5483 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5484 if (pointAngle<cpointAngle) continue;
5486 Bool_t isGamma = kFALSE;
5487 vertex.SetParamP(*track0); //track0 - plus
5488 vertex.SetParamN(*track1); //track1 - minus
5489 vertex.Update(fprimvertex);
5490 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5491 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5493 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5494 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5495 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5496 Float_t sigmae = 0.15*0.15;
5497 if (vertex.GetRr()<80)
5498 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5499 sigmae+= TMath::Sqrt(sigmae);
5500 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5501 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5502 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5503 Int_t row0 = GetRowNumber(vertex.GetRr());
5505 //Bo: if (vertex.GetDist2()>0.2) continue;
5506 if (vertex.GetDcaV0Daughters()>0.2) continue;
5507 densb0 = track0->Density2(0,row0-5);
5508 densb1 = track1->Density2(0,row0-5);
5509 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5510 densa0 = track0->Density2(row0+5,row0+40);
5511 densa1 = track1->Density2(row0+5,row0+40);
5512 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5515 densa0 = track0->Density2(0,40); //cluster density
5516 densa1 = track1->Density2(0,40); //cluster density
5517 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5519 //Bo: vertex.SetLab(0,track0->GetLabel());
5520 //Bo: vertex.SetLab(1,track1->GetLabel());
5521 vertex.SetChi2After((densa0+densa1)*0.5);
5522 vertex.SetChi2Before((densb0+densb1)*0.5);
5523 vertex.SetIndex(0,i);
5524 vertex.SetIndex(1,j);
5525 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5526 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5527 //Bo: vertex.SetRp(track0->TPCrPIDs());
5528 //Bo: vertex.SetRm(track1->TPCrPIDs());
5529 tpcv0s->AddLast(new AliESDv0(vertex));
5532 Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
5533 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5534 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5535 if (AliTPCReconstructor::StreamLevel()>1) {
5536 Int_t lab0=track0->GetLabel();
5537 Int_t lab1=track1->GetLabel();
5538 Char_t c0=track0->GetCircular();
5539 Char_t c1=track1->GetCircular();
5542 "vertex.="<<&vertex<<
5545 "Helix0.="<<&helixes[i]<<
5548 "Helix1.="<<&helixes[j]<<
5549 "pointAngle="<<pointAngle<<
5550 "pointAngle2="<<pointAngle2<<
5555 "zvertex="<<zvertex<<
5558 "npoints="<<npoints<<
5559 "radius0="<<radius[0]<<
5560 "delta0="<<delta[0]<<
5561 "radius1="<<radius[1]<<
5562 "delta1="<<delta[1]<<
5563 "radiusm="<<radiusm<<
5565 "sdcar0="<<sdcar[i]<<
5566 "sdcar1="<<sdcar[j]<<
5567 "cdcar0="<<cdcar[i]<<
5568 "cdcar1="<<cdcar[j]<<
5569 "pulldcar0="<<pulldcar[i]<<
5570 "pulldcar1="<<pulldcar[j]<<
5571 "pulldcaz0="<<pulldcaz[i]<<
5572 "pulldcaz1="<<pulldcaz[j]<<
5573 "pulldca0="<<pulldca[i]<<
5574 "pulldca1="<<pulldca[j]<<
5584 Float_t *quality = new Float_t[ncandidates];
5585 Int_t *indexes = new Int_t[ncandidates];
5587 for (Int_t i=0;i<ncandidates;i++){
5589 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5590 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5591 // quality[i] /= (0.5+v0->GetDist2());
5592 // quality[i] *= v0->GetChi2After(); //density factor
5594 Int_t index0 = v0->GetIndex(0);
5595 Int_t index1 = v0->GetIndex(1);
5596 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5597 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5601 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5602 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5603 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5604 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5607 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5610 for (Int_t i=0;i<ncandidates;i++){
5611 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5613 Int_t index0 = v0->GetIndex(0);
5614 Int_t index1 = v0->GetIndex(1);
5615 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5616 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5617 if (!track0||!track1) {
5621 Bool_t accept =kTRUE; //default accept
5622 Int_t *v0indexes0 = track0->GetV0Indexes();
5623 Int_t *v0indexes1 = track1->GetV0Indexes();
5625 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5626 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5627 if (v0indexes0[1]!=0) order0 =2;
5628 if (v0indexes1[1]!=0) order1 =2;
5630 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5631 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5633 AliESDv0 * v02 = v0;
5635 //Bo: v0->SetOrder(0,order0);
5636 //Bo: v0->SetOrder(1,order1);
5637 //Bo: v0->SetOrder(1,order0+order1);
5638 v0->SetOnFlyStatus(kTRUE);
5639 Int_t index = esd->AddV0(v0);
5640 v02 = esd->GetV0(index);
5641 v0indexes0[order0]=index;
5642 v0indexes1[order1]=index;
5646 Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
5647 if (AliTPCReconstructor::StreamLevel()>1) {
5648 Int_t lab0=track0->GetLabel();
5649 Int_t lab1=track1->GetLabel();
5658 "dca0="<<dca[index0]<<
5659 "dca1="<<dca[index1]<<
5663 "quality="<<quality[i]<<
5664 "pulldca0="<<pulldca[index0]<<
5665 "pulldca1="<<pulldca[index1]<<
5689 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5693 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5696 // refit kink towards to the vertex
5699 AliKink &kink=(AliKink &)knk;
5701 Int_t row0 = GetRowNumber(kink.GetR());
5702 FollowProlongation(mother,0);
5703 mother.Reset(kFALSE);
5705 FollowProlongation(daughter,row0);
5706 daughter.Reset(kFALSE);
5707 FollowBackProlongation(daughter,158);
5708 daughter.Reset(kFALSE);
5709 Int_t first = TMath::Max(row0-20,30);
5710 Int_t last = TMath::Min(row0+20,140);
5712 const Int_t kNdiv =5;
5713 AliTPCseed param0[kNdiv]; // parameters along the track
5714 AliTPCseed param1[kNdiv]; // parameters along the track
5715 AliKink kinks[kNdiv]; // corresponding kink parameters
5718 for (Int_t irow=0; irow<kNdiv;irow++){
5719 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5721 // store parameters along the track
5723 for (Int_t irow=0;irow<kNdiv;irow++){
5724 FollowBackProlongation(mother, rows[irow]);
5725 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5726 param0[irow] = mother;
5727 param1[kNdiv-1-irow] = daughter;
5731 for (Int_t irow=0; irow<kNdiv-1;irow++){
5732 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5733 kinks[irow].SetMother(param0[irow]);
5734 kinks[irow].SetDaughter(param1[irow]);
5735 kinks[irow].Update();
5738 // choose kink with best "quality"
5740 Double_t mindist = 10000;
5741 for (Int_t irow=0;irow<kNdiv;irow++){
5742 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5743 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5744 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5746 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5747 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5748 if (normdist < mindist){
5754 if (index==-1) return 0;
5757 param0[index].Reset(kTRUE);
5758 FollowProlongation(param0[index],0);
5760 mother = param0[index];
5761 daughter = param1[index]; // daughter in vertex
5763 kink.SetMother(mother);
5764 kink.SetDaughter(daughter);
5766 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5767 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5768 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5769 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5770 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5771 mother.SetLabel(kink.GetLabel(0));
5772 daughter.SetLabel(kink.GetLabel(1));
5778 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5780 // update Kink quality information for mother after back propagation
5782 if (seed->GetKinkIndex(0)>=0) return;
5783 for (Int_t ikink=0;ikink<3;ikink++){
5784 Int_t index = seed->GetKinkIndex(ikink);
5785 if (index>=0) break;
5786 index = TMath::Abs(index)-1;
5787 AliESDkink * kink = fEvent->GetKink(index);
5788 kink->SetTPCDensity(-1,0,0);
5789 kink->SetTPCDensity(1,0,1);
5791 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5792 if (row0<15) row0=15;
5794 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5795 if (row1>145) row1=145;
5797 Int_t found,foundable,shared;
5798 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5799 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5800 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5801 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5806 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5808 // update Kink quality information for daughter after refit
5810 if (seed->GetKinkIndex(0)<=0) return;
5811 for (Int_t ikink=0;ikink<3;ikink++){
5812 Int_t index = seed->GetKinkIndex(ikink);
5813 if (index<=0) break;
5814 index = TMath::Abs(index)-1;
5815 AliESDkink * kink = fEvent->GetKink(index);
5816 kink->SetTPCDensity(-1,1,0);
5817 kink->SetTPCDensity(-1,1,1);
5819 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5820 if (row0<15) row0=15;
5822 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5823 if (row1>145) row1=145;
5825 Int_t found,foundable,shared;
5826 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5827 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5828 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5829 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5835 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5838 // check kink point for given track
5839 // if return value=0 kink point not found
5840 // otherwise seed0 correspond to mother particle
5841 // seed1 correspond to daughter particle
5842 // kink parameter of kink point
5843 AliKink &kink=(AliKink &)knk;
5845 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5846 Int_t first = seed->GetFirstPoint();
5847 Int_t last = seed->GetLastPoint();
5848 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5851 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5852 if (!seed1) return 0;
5853 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5854 seed1->Reset(kTRUE);
5855 FollowProlongation(*seed1,158);
5856 seed1->Reset(kTRUE);
5857 last = seed1->GetLastPoint();
5859 AliTPCseed *seed0 = new AliTPCseed(*seed);
5860 seed0->Reset(kFALSE);
5863 AliTPCseed param0[20]; // parameters along the track
5864 AliTPCseed param1[20]; // parameters along the track
5865 AliKink kinks[20]; // corresponding kink parameters
5867 for (Int_t irow=0; irow<20;irow++){
5868 rows[irow] = first +((last-first)*irow)/19;
5870 // store parameters along the track
5872 for (Int_t irow=0;irow<20;irow++){
5873 FollowBackProlongation(*seed0, rows[irow]);
5874 FollowProlongation(*seed1,rows[19-irow]);
5875 param0[irow] = *seed0;
5876 param1[19-irow] = *seed1;
5880 for (Int_t irow=0; irow<19;irow++){
5881 kinks[irow].SetMother(param0[irow]);
5882 kinks[irow].SetDaughter(param1[irow]);
5883 kinks[irow].Update();
5886 // choose kink with biggest change of angle
5888 Double_t maxchange= 0;
5889 for (Int_t irow=1;irow<19;irow++){
5890 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5891 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5892 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5893 if ( quality > maxchange){
5894 maxchange = quality;
5901 if (index<0) return 0;
5903 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5904 seed0 = new AliTPCseed(param0[index]);
5905 seed1 = new AliTPCseed(param1[index]);
5906 seed0->Reset(kFALSE);
5907 seed1->Reset(kFALSE);
5908 seed0->ResetCovariance(10.);
5909 seed1->ResetCovariance(10.);
5910 FollowProlongation(*seed0,0);
5911 FollowBackProlongation(*seed1,158);
5912 mother = *seed0; // backup mother at position 0
5913 seed0->Reset(kFALSE);
5914 seed1->Reset(kFALSE);
5915 seed0->ResetCovariance(10.);
5916 seed1->ResetCovariance(10.);
5918 first = TMath::Max(row0-20,0);
5919 last = TMath::Min(row0+20,158);
5921 for (Int_t irow=0; irow<20;irow++){
5922 rows[irow] = first +((last-first)*irow)/19;
5924 // store parameters along the track
5926 for (Int_t irow=0;irow<20;irow++){
5927 FollowBackProlongation(*seed0, rows[irow]);
5928 FollowProlongation(*seed1,rows[19-irow]);
5929 param0[irow] = *seed0;
5930 param1[19-irow] = *seed1;
5934 for (Int_t irow=0; irow<19;irow++){
5935 kinks[irow].SetMother(param0[irow]);
5936 kinks[irow].SetDaughter(param1[irow]);
5937 // param0[irow].Dump();
5938 //param1[irow].Dump();
5939 kinks[irow].Update();
5942 // choose kink with biggest change of angle
5945 for (Int_t irow=0;irow<20;irow++){
5946 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5947 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5948 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5949 if ( quality > maxchange){
5950 maxchange = quality;
5957 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5962 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5964 kink.SetMother(param0[index]);
5965 kink.SetDaughter(param1[index]);
5967 row0 = GetRowNumber(kink.GetR());
5968 kink.SetTPCRow0(row0);
5969 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5970 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5971 kink.SetIndex(-10,0);
5972 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5973 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5974 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5977 // new (&mother) AliTPCseed(param0[index]);
5978 daughter = param1[index];
5979 daughter.SetLabel(kink.GetLabel(1));
5980 param0[index].Reset(kTRUE);
5981 FollowProlongation(param0[index],0);
5982 mother = param0[index];
5983 mother.SetLabel(kink.GetLabel(0));
5993 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5996 // reseed - refit - track
5999 // Int_t last = fSectors->GetNRows()-1;
6001 if (fSectors == fOuterSec){
6002 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6006 first = t->GetFirstPoint();
6008 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6009 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6011 FollowProlongation(*t,first);
6021 //_____________________________________________________________________________
6022 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6023 //-----------------------------------------------------------------
6024 // This function reades track seeds.
6025 //-----------------------------------------------------------------
6026 TDirectory *savedir=gDirectory;
6028 TFile *in=(TFile*)inp;
6029 if (!in->IsOpen()) {
6030 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6035 TTree *seedTree=(TTree*)in->Get("Seeds");
6037 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6038 cerr<<"can't get a tree with track seeds !\n";
6041 AliTPCtrack *seed=new AliTPCtrack;
6042 seedTree->SetBranchAddress("tracks",&seed);
6044 if (fSeeds==0) fSeeds=new TObjArray(15000);
6046 Int_t n=(Int_t)seedTree->GetEntries();
6047 for (Int_t i=0; i<n; i++) {
6048 seedTree->GetEvent(i);
6049 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6058 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6061 if (fSeeds) DeleteSeeds();
6064 if (!fSeeds) return 1;
6071 //_____________________________________________________________________________
6072 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6073 //-----------------------------------------------------------------
6074 // This is a track finder.
6075 //-----------------------------------------------------------------
6076 TDirectory *savedir=gDirectory;
6080 fSeeds = Tracking();
6083 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6085 //activate again some tracks
6086 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6087 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6089 Int_t nc=t.GetNumberOfClusters();
6091 delete fSeeds->RemoveAt(i);
6095 if (pt->GetRemoval()==10) {
6096 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6097 pt->Desactivate(10); // make track again active
6099 pt->Desactivate(20);
6100 delete fSeeds->RemoveAt(i);
6105 RemoveUsed2(fSeeds,0.85,0.85,0);
6106 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6107 //FindCurling(fSeeds, fEvent,0);
6108 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6109 RemoveUsed2(fSeeds,0.5,0.4,20);
6110 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6113 // // refit short tracks
6115 Int_t nseed=fSeeds->GetEntriesFast();
6118 for (Int_t i=0; i<nseed; i++) {
6119 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6121 Int_t nc=t.GetNumberOfClusters();
6123 delete fSeeds->RemoveAt(i);
6126 CookLabel(pt,0.1); //For comparison only
6127 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6128 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6130 if (fDebug>0) cerr<<found<<'\r';
6134 delete fSeeds->RemoveAt(i);
6138 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6140 //RemoveUsed(fSeeds,0.9,0.9,6);
6142 nseed=fSeeds->GetEntriesFast();
6144 for (Int_t i=0; i<nseed; i++) {
6145 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6147 Int_t nc=t.GetNumberOfClusters();
6149 delete fSeeds->RemoveAt(i);
6153 t.CookdEdx(0.02,0.6);
6154 // CheckKinkPoint(&t,0.05);
6155 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6156 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6164 delete fSeeds->RemoveAt(i);
6165 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6167 // FollowProlongation(*seed1,0);
6168 // Int_t n = seed1->GetNumberOfClusters();
6169 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6170 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6173 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6177 SortTracks(fSeeds, 1);
6181 PrepareForBackProlongation(fSeeds,5.);
6182 PropagateBack(fSeeds);
6183 printf("Time for back propagation: \t");timer.Print();timer.Start();
6187 PrepareForProlongation(fSeeds,5.);
6188 PropagateForward2(fSeeds);
6190 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6191 // RemoveUsed(fSeeds,0.7,0.7,6);
6192 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6194 nseed=fSeeds->GetEntriesFast();
6196 for (Int_t i=0; i<nseed; i++) {
6197 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6199 Int_t nc=t.GetNumberOfClusters();
6201 delete fSeeds->RemoveAt(i);
6204 t.CookdEdx(0.02,0.6);
6205 // CookLabel(pt,0.1); //For comparison only
6206 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6207 if ((pt->IsActive() || (pt->fRemoval==10) )){
6208 cerr<<found++<<'\r';
6211 delete fSeeds->RemoveAt(i);
6216 // fNTracks = found;
6218 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6221 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6222 Info("Clusters2Tracks","Number of found tracks %d",found);
6224 // UnloadClusters();
6229 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6232 // tracking of the seeds
6235 fSectors = fOuterSec;
6236 ParallelTracking(arr,150,63);
6237 fSectors = fOuterSec;
6238 ParallelTracking(arr,63,0);
6241 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6246 TObjArray * arr = new TObjArray;
6248 fSectors = fOuterSec;
6251 for (Int_t sec=0;sec<fkNOS;sec++){
6252 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6253 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6254 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6257 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6269 TObjArray * AliTPCtrackerMI::Tracking()
6273 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6276 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6278 TObjArray * seeds = new TObjArray;
6287 Float_t fnumber = 3.0;
6288 Float_t fdensity = 3.0;
6293 for (Int_t delta = 0; delta<18; delta+=6){
6297 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6298 SumTracks(seeds,arr);
6299 SignClusters(seeds,fnumber,fdensity);
6301 for (Int_t i=2;i<6;i+=2){
6302 // seed high pt tracks
6305 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6306 SumTracks(seeds,arr);
6307 SignClusters(seeds,fnumber,fdensity);
6312 // RemoveUsed(seeds,0.9,0.9,1);
6313 // UnsignClusters();
6314 // SignClusters(seeds,fnumber,fdensity);
6318 for (Int_t delta = 20; delta<120; delta+=10){
6320 // seed high pt tracks
6324 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6325 SumTracks(seeds,arr);
6326 SignClusters(seeds,fnumber,fdensity);
6331 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6332 SumTracks(seeds,arr);
6333 SignClusters(seeds,fnumber,fdensity);
6344 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6348 // RemoveUsed(seeds,0.75,0.75,1);
6350 //SignClusters(seeds,fnumber,fdensity);
6359 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6360 SumTracks(seeds,arr);
6361 SignClusters(seeds,fnumber,fdensity);
6363 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6364 SumTracks(seeds,arr);
6365 SignClusters(seeds,fnumber,fdensity);
6367 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6368 SumTracks(seeds,arr);
6369 SignClusters(seeds,fnumber,fdensity);
6373 for (Int_t delta = 3; delta<30; delta+=5){
6379 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6380 SumTracks(seeds,arr);
6381 SignClusters(seeds,fnumber,fdensity);
6383 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6384 SumTracks(seeds,arr);
6385 SignClusters(seeds,fnumber,fdensity);
6396 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6399 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6405 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6406 SumTracks(seeds,arr);
6407 SignClusters(seeds,fnumber,fdensity);
6409 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6410 SumTracks(seeds,arr);
6411 SignClusters(seeds,fnumber,fdensity);
6415 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6426 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6429 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6430 // no primary vertex seeding tried
6434 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6436 TObjArray * seeds = new TObjArray;
6441 Float_t fnumber = 3.0;
6442 Float_t fdensity = 3.0;
6445 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6446 cuts[1] = 3.5; // max tan(phi) angle for seeding
6447 cuts[2] = 3.; // not used (cut on z primary vertex)
6448 cuts[3] = 3.5; // max tan(theta) angle for seeding
6450 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6452 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6453 SumTracks(seeds,arr);
6454 SignClusters(seeds,fnumber,fdensity);
6458 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6469 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6472 //sum tracks to common container
6473 //remove suspicious tracks
6474 Int_t nseed = arr2->GetEntriesFast();
6475 for (Int_t i=0;i<nseed;i++){
6476 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6479 // remove tracks with too big curvature
6481 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6482 delete arr2->RemoveAt(i);
6485 // REMOVE VERY SHORT TRACKS
6486 if (pt->GetNumberOfClusters()<20){
6487 delete arr2->RemoveAt(i);
6490 // NORMAL ACTIVE TRACK
6491 if (pt->IsActive()){
6492 arr1->AddLast(arr2->RemoveAt(i));
6495 //remove not usable tracks
6496 if (pt->GetRemoval()!=10){
6497 delete arr2->RemoveAt(i);
6501 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6502 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6503 arr1->AddLast(arr2->RemoveAt(i));
6505 delete arr2->RemoveAt(i);
6509 delete arr2; arr2 = 0;
6514 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6517 // try to track in parralel
6519 Int_t nseed=arr->GetEntriesFast();
6520 //prepare seeds for tracking
6521 for (Int_t i=0; i<nseed; i++) {
6522 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6524 if (!t.IsActive()) continue;
6525 // follow prolongation to the first layer
6526 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6527 FollowProlongation(t, rfirst+1);
6532 for (Int_t nr=rfirst; nr>=rlast; nr--){
6533 if (nr<fInnerSec->GetNRows())
6534 fSectors = fInnerSec;
6536 fSectors = fOuterSec;
6537 // make indexes with the cluster tracks for given
6539 // find nearest cluster
6540 for (Int_t i=0; i<nseed; i++) {
6541 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6543 if (nr==80) pt->UpdateReference();
6544 if (!pt->IsActive()) continue;
6545 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6546 if (pt->GetRelativeSector()>17) {
6549 UpdateClusters(t,nr);
6551 // prolonagate to the nearest cluster - if founded
6552 for (Int_t i=0; i<nseed; i++) {
6553 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6555 if (!pt->IsActive()) continue;
6556 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6557 if (pt->GetRelativeSector()>17) {
6560 FollowToNextCluster(*pt,nr);
6565 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6569 // if we use TPC track itself we have to "update" covariance
6571 Int_t nseed= arr->GetEntriesFast();
6572 for (Int_t i=0;i<nseed;i++){
6573 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6577 //rotate to current local system at first accepted point
6578 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6579 Int_t sec = (index&0xff000000)>>24;
6581 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6582 if (angle1>TMath::Pi())
6583 angle1-=2.*TMath::Pi();
6584 Float_t angle2 = pt->GetAlpha();
6586 if (TMath::Abs(angle1-angle2)>0.001){
6587 pt->Rotate(angle1-angle2);
6588 //angle2 = pt->GetAlpha();
6589 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6590 //if (pt->GetAlpha()<0)
6591 // pt->fRelativeSector+=18;
6592 //sec = pt->fRelativeSector;
6601 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6605 // if we use TPC track itself we have to "update" covariance
6607 Int_t nseed= arr->GetEntriesFast();
6608 for (Int_t i=0;i<nseed;i++){
6609 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6612 pt->SetFirstPoint(pt->GetLastPoint());
6620 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6623 // make back propagation
6625 Int_t nseed= arr->GetEntriesFast();
6626 for (Int_t i=0;i<nseed;i++){
6627 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6628 if (pt&& pt->GetKinkIndex(0)<=0) {
6629 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6630 fSectors = fInnerSec;
6631 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6632 //fSectors = fOuterSec;
6633 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6634 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6635 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6636 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6639 if (pt&& pt->GetKinkIndex(0)>0) {
6640 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6641 pt->SetFirstPoint(kink->GetTPCRow0());
6642 fSectors = fInnerSec;
6643 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6651 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6654 // make forward propagation
6656 Int_t nseed= arr->GetEntriesFast();
6658 for (Int_t i=0;i<nseed;i++){
6659 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6661 FollowProlongation(*pt,0);
6670 Int_t AliTPCtrackerMI::PropagateForward()
6673 // propagate track forward
6675 Int_t nseed = fSeeds->GetEntriesFast();
6676 for (Int_t i=0;i<nseed;i++){
6677 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6679 AliTPCseed &t = *pt;
6680 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6681 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6682 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6683 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6687 fSectors = fOuterSec;
6688 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6689 fSectors = fInnerSec;
6690 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6700 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6703 // make back propagation, in between row0 and row1
6707 fSectors = fInnerSec;
6710 if (row1<fSectors->GetNRows())
6713 r1 = fSectors->GetNRows()-1;
6715 if (row0<fSectors->GetNRows()&& r1>0 )
6716 FollowBackProlongation(*pt,r1);
6717 if (row1<=fSectors->GetNRows())
6720 r1 = row1 - fSectors->GetNRows();
6721 if (r1<=0) return 0;
6722 if (r1>=fOuterSec->GetNRows()) return 0;
6723 fSectors = fOuterSec;
6724 return FollowBackProlongation(*pt,r1);
6732 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6736 Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6737 // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6738 Float_t padlength = GetPadPitchLength(row);
6740 Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6741 Double_t angulary = seed->GetSnp();
6742 angulary = angulary*angulary/(1-angulary*angulary);
6743 seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6745 Float_t sresz = fParam->GetZSigma();
6746 Float_t angularz = seed->GetTgl();
6747 seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6749 Float_t wy = GetSigmaY(seed);
6750 Float_t wz = GetSigmaZ(seed);
6753 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6754 printf("problem\n");
6760 Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
6764 Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6765 Float_t padlength = fParam->GetPadPitchLength(seed->GetSector());
6766 Float_t sres = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6767 Float_t angular = seed->GetSnp();
6768 angular = angular*angular/(1-angular*angular);
6769 // angular*=angular;
6770 //angular = TMath::Sqrt(angular/(1-angular));
6771 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
6774 Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
6778 Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6779 Float_t padlength = fParam->GetPadPitchLength(seed->GetSector());
6780 Float_t sres = fParam->GetZSigma();
6781 Float_t angular = seed->GetTgl();
6782 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
6787 //__________________________________________________________________________
6788 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6789 //--------------------------------------------------------------------
6790 //This function "cooks" a track label. If label<0, this track is fake.
6791 //--------------------------------------------------------------------
6792 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6794 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6798 Int_t noc=t->GetNumberOfClusters();
6800 //printf("\nnot founded prolongation\n\n\n");
6806 AliTPCclusterMI *clusters[160];
6808 for (Int_t i=0;i<160;i++) {
6815 for (i=0; i<160 && current<noc; i++) {
6817 Int_t index=t->GetClusterIndex2(i);
6818 if (index<=0) continue;
6819 if (index&0x8000) continue;
6821 //clusters[current]=GetClusterMI(index);
6822 if (t->GetClusterPointer(i)){
6823 clusters[current]=t->GetClusterPointer(i);
6829 Int_t lab=123456789;
6830 for (i=0; i<noc; i++) {
6831 AliTPCclusterMI *c=clusters[i];
6833 lab=TMath::Abs(c->GetLabel(0));
6835 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6841 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6843 for (i=0; i<noc; i++) {
6844 AliTPCclusterMI *c=clusters[i];
6846 if (TMath::Abs(c->GetLabel(1)) == lab ||
6847 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6850 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6853 Int_t tail=Int_t(0.10*noc);
6856 for (i=1; i<=160&&ind<tail; i++) {
6857 // AliTPCclusterMI *c=clusters[noc-i];
6858 AliTPCclusterMI *c=clusters[i];
6860 if (lab == TMath::Abs(c->GetLabel(0)) ||
6861 lab == TMath::Abs(c->GetLabel(1)) ||
6862 lab == TMath::Abs(c->GetLabel(2))) max++;
6865 if (max < Int_t(0.5*tail)) lab=-lab;
6872 //delete[] clusters;
6876 //__________________________________________________________________________
6877 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6878 //--------------------------------------------------------------------
6879 //This function "cooks" a track label. If label<0, this track is fake.
6880 //--------------------------------------------------------------------
6881 Int_t noc=t->GetNumberOfClusters();
6883 //printf("\nnot founded prolongation\n\n\n");
6889 AliTPCclusterMI *clusters[160];
6891 for (Int_t i=0;i<160;i++) {
6898 for (i=0; i<160 && current<noc; i++) {
6899 if (i<first) continue;
6900 if (i>last) continue;
6901 Int_t index=t->GetClusterIndex2(i);
6902 if (index<=0) continue;
6903 if (index&0x8000) continue;
6905 //clusters[current]=GetClusterMI(index);
6906 if (t->GetClusterPointer(i)){
6907 clusters[current]=t->GetClusterPointer(i);
6912 if (noc<5) return -1;
6913 Int_t lab=123456789;
6914 for (i=0; i<noc; i++) {
6915 AliTPCclusterMI *c=clusters[i];
6917 lab=TMath::Abs(c->GetLabel(0));
6919 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6925 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6927 for (i=0; i<noc; i++) {
6928 AliTPCclusterMI *c=clusters[i];
6930 if (TMath::Abs(c->GetLabel(1)) == lab ||
6931 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6934 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6937 Int_t tail=Int_t(0.10*noc);
6940 for (i=1; i<=160&&ind<tail; i++) {
6941 // AliTPCclusterMI *c=clusters[noc-i];
6942 AliTPCclusterMI *c=clusters[i];
6944 if (lab == TMath::Abs(c->GetLabel(0)) ||
6945 lab == TMath::Abs(c->GetLabel(1)) ||
6946 lab == TMath::Abs(c->GetLabel(2))) max++;
6949 if (max < Int_t(0.5*tail)) lab=-lab;
6952 // t->SetLabel(lab);
6956 //delete[] clusters;
6960 Int_t AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const
6962 //return pad row number for this x
6965 r=fRow[fN-1].GetX();
6966 if (x > r) return fN;
6968 if (x < r) return -1;
6969 return Int_t((x-r)/fPadPitchLength + 0.5);}
6971 r=fRow[fN-1].GetX();
6972 if (x > r) return fN;
6974 if (x < r) return -1;
6975 Double_t r1=fRow[64].GetX();
6977 return Int_t((x-r)/f1PadPitchLength + 0.5);}
6979 return (Int_t((x-r1)/f2PadPitchLength + 0.5)+64);}
6983 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6985 //return pad row number for given x vector
6986 Float_t phi = TMath::ATan2(x[1],x[0]);
6987 if(phi<0) phi=2.*TMath::Pi()+phi;
6988 // Get the local angle in the sector philoc
6989 const Float_t kRaddeg = 180/3.14159265358979312;
6990 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6991 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6992 return GetRowNumber(localx);
6995 //_________________________________________________________________________
6996 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
6997 //-----------------------------------------------------------------------
6998 // Setup inner sector
6999 //-----------------------------------------------------------------------
7001 fAlpha=par->GetInnerAngle();
7002 fAlphaShift=par->GetInnerAngleShift();
7003 fPadPitchWidth=par->GetInnerPadPitchWidth();
7004 fPadPitchLength=par->GetInnerPadPitchLength();
7005 fN=par->GetNRowLow();
7006 if(fRow)delete [] fRow;fRow = 0;
7007 fRow=new AliTPCRow[fN];
7008 for (Int_t i=0; i<fN; i++) {
7009 fRow[i].SetX(par->GetPadRowRadiiLow(i));
7010 fRow[i].SetDeadZone(1.5); //1.5 cm of dead zone
7013 fAlpha=par->GetOuterAngle();
7014 fAlphaShift=par->GetOuterAngleShift();
7015 fPadPitchWidth = par->GetOuterPadPitchWidth();
7016 fPadPitchLength = par->GetOuter1PadPitchLength();
7017 f1PadPitchLength = par->GetOuter1PadPitchLength();
7018 f2PadPitchLength = par->GetOuter2PadPitchLength();
7019 fN=par->GetNRowUp();
7020 if(fRow)delete [] fRow;fRow = 0;
7021 fRow=new AliTPCRow[fN];
7022 for (Int_t i=0; i<fN; i++) {
7023 fRow[i].SetX(par->GetPadRowRadiiUp(i));
7024 fRow[i].SetDeadZone(1.5); // 1.5 cm of dead zone
7029 AliTPCtrackerMI::AliTPCRow::AliTPCRow():
7039 // default constructor
7043 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
7045 for (Int_t i = 0; i < fN1; i++)
7046 fClusters1[i].~AliTPCclusterMI();
7047 delete [] fClusters1;
7048 for (Int_t i = 0; i < fN2; i++)
7049 fClusters2[i].~AliTPCclusterMI();
7050 delete [] fClusters2;
7055 //_________________________________________________________________________
7057 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
7058 //-----------------------------------------------------------------------
7059 // Insert a cluster into this pad row in accordence with its y-coordinate
7060 //-----------------------------------------------------------------------
7061 if (fN==kMaxClusterPerRow) {
7062 //AliInfo("AliTPCRow::InsertCluster(): Too many clusters");
7066 //AliInfo("AliTPCRow::InsertCluster(): Too many clusters !");
7069 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
7070 Int_t i=Find(c->GetZ());
7071 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
7072 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
7073 fIndex[i]=index; fClusters[i]=c; fN++;
7076 void AliTPCtrackerMI::AliTPCRow::ResetClusters() {
7079 // MvL: Need to call destructors for AliTPCclusterMI, to delete fInfo
7080 for (Int_t i = 0; i < fN1; i++)
7081 fClusters1[i].~AliTPCclusterMI();
7082 delete [] fClusters1; fClusters1=0;
7083 for (Int_t i = 0; i < fN2; i++)
7084 fClusters2[i].~AliTPCclusterMI();
7085 delete [] fClusters2; fClusters2=0;
7090 //delete[] fClusterArray;
7096 //___________________________________________________________________
7097 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
7098 //-----------------------------------------------------------------------
7099 // Return the index of the nearest cluster
7100 //-----------------------------------------------------------------------
7101 if (fN==0) return 0;
7102 if (z <= fClusters[0]->GetZ()) return 0;
7103 if (z > fClusters[fN-1]->GetZ()) return fN;
7104 Int_t b=0, e=fN-1, m=(b+e)/2;
7105 for (; b<e; m=(b+e)/2) {
7106 if (z > fClusters[m]->GetZ()) b=m+1;
7114 //___________________________________________________________________
7115 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
7116 //-----------------------------------------------------------------------
7117 // Return the index of the nearest cluster in z y
7118 //-----------------------------------------------------------------------
7119 Float_t maxdistance = roady*roady + roadz*roadz;
7121 AliTPCclusterMI *cl =0;
7122 for (Int_t i=Find(z-roadz); i<fN; i++) {
7123 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
7124 if (c->GetZ() > z+roadz) break;
7125 if ( (c->GetY()-y) > roady ) continue;
7126 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
7127 if (maxdistance>distance) {
7128 maxdistance = distance;
7135 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
7137 //-----------------------------------------------------------------------
7138 // Return the index of the nearest cluster in z y
7139 //-----------------------------------------------------------------------
7140 Float_t maxdistance = roady*roady + roadz*roadz;
7141 AliTPCclusterMI *cl =0;
7143 //PH Check boundaries. 510 is the size of fFastCluster
7144 Int_t iz1 = Int_t(z-roadz+254.5);
7145 if (iz1<0 || iz1>=510) return cl;
7146 iz1 = TMath::Max(GetFastCluster(iz1)-1,0);
7147 Int_t iz2 = Int_t(z+roadz+255.5);
7148 if (iz2<0 || iz2>=510) return cl;
7149 iz2 = TMath::Min(GetFastCluster(iz2)+1,fN);
7150 Bool_t skipUsed = !(AliTPCReconstructor::GetRecoParam()->GetClusterSharing());
7151 //FindNearest3(y,z,roady,roadz,index);
7152 // for (Int_t i=Find(z-roadz); i<fN; i++) {
7153 for (Int_t i=iz1; i<iz2; i++) {
7154 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
7155 if (c->GetZ() > z+roadz) break;
7156 if ( c->GetY()-y > roady ) continue;
7157 if ( y-c->GetY() > roady ) continue;
7158 if (skipUsed && c->IsUsed(11)) continue;
7159 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
7160 if (maxdistance>distance) {
7161 maxdistance = distance;
7164 //roady = TMath::Sqrt(maxdistance);
7171 void AliTPCtrackerMI::AliTPCRow::SetFastCluster(Int_t i, Short_t cl){
7173 // Set cluster info for fast navigation
7184 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
7186 //-----------------------------------------------------------------------
7187 // Fill the cluster and sharing bitmaps of the track
7188 //-----------------------------------------------------------------------
7190 Int_t firstpoint = 0;
7191 Int_t lastpoint = 159;
7192 AliTPCTrackerPoint *point;
7194 for (int iter=firstpoint; iter<lastpoint; iter++) {
7195 point = t->GetTrackPoint(iter);
7197 t->SetClusterMapBit(iter, kTRUE);
7198 if (point->IsShared())
7199 t->SetSharedMapBit(iter,kTRUE);
7201 t->SetSharedMapBit(iter, kFALSE);
7204 t->SetClusterMapBit(iter, kFALSE);
7205 t->SetSharedMapBit(iter, kFALSE);