Removing error message in SetFastCluster (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerMI.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 //-------------------------------------------------------
18 //          Implementation of the TPC tracker
19 //
20 //   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
21 // 
22 //  AliTPC parallel tracker
23 //
24 //  The track fitting is based on Kalaman filtering approach
25
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
31
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) ...
35
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) 
39
40 //  Repeating several seeding - tracking procedures some of the tracks can be find 
41 //  several times. 
42
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
53 //  
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 
59 //                  low efficiency
60 //                 
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)
66 //                          
67 //
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
71 //
72 // The debug level -  different procedure produce tree for numerical debugging
73 //                    To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
74 //
75
76 //-------------------------------------------------------
77
78
79 /* $Id$ */
80
81 #include "Riostream.h"
82 #include <TClonesArray.h>
83 #include <TFile.h>
84 #include <TObjArray.h>
85 #include <TTree.h>
86 #include "AliLog.h"
87 #include "AliComplexCluster.h"
88 #include "AliESDEvent.h"
89 #include "AliESDtrack.h"
90 #include "AliESDVertex.h"
91 #include "AliKink.h"
92 #include "AliV0.h"
93 #include "AliHelix.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"
104 #include "AliPID.h"
105 #include "TTreeStream.h"
106 #include "AliAlignObj.h"
107 #include "AliTrackPointArray.h"
108 #include "TRandom.h"
109 #include "AliTPCcalibDB.h"
110 #include "AliTPCTransform.h"
111
112 //
113
114 ClassImp(AliTPCtrackerMI)
115
116
117 class AliTPCFastMath {
118 public:
119   AliTPCFastMath();  
120   static Double_t FastAsin(Double_t x);   
121  private: 
122   static Double_t fgFastAsin[20000];  //lookup table for fast asin computation
123 };
124
125 Double_t AliTPCFastMath::fgFastAsin[20000];
126 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
127
128 AliTPCFastMath::AliTPCFastMath(){
129   //
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]);
134   }
135 }
136
137 Double_t AliTPCFastMath::FastAsin(Double_t x){
138   //
139   // return asin using lookup table
140   if (x>0){
141     Int_t index = int(x*10000);
142     return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
143   }
144   x*=-1;
145   Int_t index = int(x*10000);
146   return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
147 }
148 //__________________________________________________________________
149 AliTPCtrackerMI::AliTPCtrackerMI()
150                 :AliTracker(),
151                  fkNIS(0),
152                  fInnerSec(0),
153                  fkNOS(0),
154                  fOuterSec(0),
155                  fN(0),
156                  fSectors(0),
157                  fInput(0),
158                  fOutput(0),
159                  fSeedTree(0),
160                  fTreeDebug(0),
161                  fEvent(0),
162                  fDebug(0),
163                  fNewIO(kFALSE),
164                  fNtracks(0),
165                  fSeeds(0),
166                  fIteration(0),
167                  fParam(0),
168                  fDebugStreamer(0)
169 {
170   //
171   // default constructor
172   //
173 }
174 //_____________________________________________________________________
175
176
177
178 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
179   //
180   //update track information using current cluster - track->fCurrentCluster
181
182
183   AliTPCclusterMI* c =track->GetCurrentCluster();
184   if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);  //sign not accepted clusters
185
186   UInt_t i = track->GetCurrentClusterIndex1();
187
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");
199   //}
200   if (track->GetFirstPoint()>track->GetRow()) 
201     track->SetFirstPoint(track->GetRow());
202   if (track->GetLastPoint()<track->GetRow()) 
203     track->SetLastPoint(track->GetRow());
204   
205
206   track->SetClusterPointer(track->GetRow(),c);  
207   //
208
209   Double_t angle2 = track->GetSnp()*track->GetSnp();
210   //
211   //SET NEW Track Point
212   //
213   if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
214   {
215     angle2 = TMath::Sqrt(angle2/(1-angle2)); 
216     AliTPCTrackerPoint   &point =*(track->GetTrackPoint(track->GetRow()));
217     //
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()));
222     //
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);
231     }
232   }  
233
234   Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
235   //
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);      
240     //}
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) ;    
247     //}
248   }
249   track->SetNoCluster(0);
250   return track->Update(c,chi2,i);
251 }
252
253
254
255 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor, 
256                                       Float_t cory, Float_t corz)
257 {
258   //
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;
265   
266   Double_t sdistancey2 = sy2+seed->GetSigmaY2();
267   Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
268   
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;
273   
274   Double_t rdistance2  = rdistancey2+rdistancez2;
275   //Int_t  accept =0;
276   
277   if (rdistance2>16) return 3;
278   
279   
280   if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)  
281     return 2;  //suspisiouce - will be changed
282   
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
286   
287   if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor ) 
288        && cluster->GetType()<0){
289     seed->SetNFoundable(seed->GetNFoundable()-1);
290     return 2;    
291   }
292   return 0;
293 }
294
295
296
297
298 //_____________________________________________________________________________
299 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par): 
300 AliTracker(), 
301                  fkNIS(par->GetNInnerSector()/2),
302                  fInnerSec(0),
303                  fkNOS(par->GetNOuterSector()/2),
304                  fOuterSec(0),
305                  fN(0),
306                  fSectors(0),
307                  fInput(0),
308                  fOutput(0),
309                  fSeedTree(0),
310                  fTreeDebug(0),
311                  fEvent(0),
312                  fDebug(0),
313                  fNewIO(0),
314                  fNtracks(0),
315                  fSeeds(0),
316                  fIteration(0),
317                  fParam(0),
318                  fDebugStreamer(0)
319 {
320   //---------------------------------------------------------------------
321   // The main TPC tracker constructor
322   //---------------------------------------------------------------------
323   fInnerSec=new AliTPCSector[fkNIS];         
324   fOuterSec=new AliTPCSector[fkNOS];
325  
326   Int_t i;
327   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
328   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
329
330   fParam = par;  
331   Int_t nrowlow = par->GetNRowLow();
332   Int_t nrowup = par->GetNRowUp();
333
334   
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());
339   }
340
341   
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());
346   }
347
348   fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
349 }
350 //________________________________________________________________________
351 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
352   AliTracker(t),
353                  fkNIS(t.fkNIS),
354                  fInnerSec(0),
355                  fkNOS(t.fkNOS),
356                  fOuterSec(0),
357                  fN(0),
358                  fSectors(0),
359                  fInput(0),
360                  fOutput(0),
361                  fSeedTree(0),
362                  fTreeDebug(0),
363                  fEvent(0),
364                  fDebug(0),
365                  fNewIO(kFALSE),
366                  fNtracks(0),
367                  fSeeds(0),
368                  fIteration(0),
369                  fParam(0),
370                  fDebugStreamer(0)
371 {
372   //------------------------------------
373   // dummy copy constructor
374   //------------------------------------------------------------------
375   fOutput=t.fOutput;
376 }
377 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
378   //------------------------------
379   // dummy 
380   //--------------------------------------------------------------
381   return *this;
382 }
383 //_____________________________________________________________________________
384 AliTPCtrackerMI::~AliTPCtrackerMI() {
385   //------------------------------------------------------------------
386   // TPC tracker destructor
387   //------------------------------------------------------------------
388   delete[] fInnerSec;
389   delete[] fOuterSec;
390   if (fSeeds) {
391     fSeeds->Delete(); 
392     delete fSeeds;
393   }
394   if (fDebugStreamer) delete fDebugStreamer;
395 }
396
397 void AliTPCtrackerMI::SetIO()
398 {
399   //
400   fNewIO   =  kTRUE;
401   fInput   =  AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName());
402   
403   fOutput  =  AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName());
404   if (fOutput){
405     AliTPCtrack *iotrack= new AliTPCtrack;
406     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
407     delete iotrack;
408   }
409 }
410
411
412 void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESDEvent * event)
413 {
414
415   // set input
416   fNewIO = kFALSE;
417   fInput    = 0;
418   fOutput   = 0;
419   fSeedTree = 0;
420   fTreeDebug =0;
421   fInput = input;
422   if (input==0){
423     return;
424   }  
425   //set output
426   fOutput = output;
427   if (output){
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);
432     delete iotrack;
433   }
434   if (output && (fDebug&2)){
435     //write the full seed information if specified in debug mode
436     //
437     fSeedTree =  new TTree("Seeds","Seeds");
438     AliTPCseed * vseed = new AliTPCseed;
439     //
440     TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
441     arrtr->ExpandCreateFast(160);
442     TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
443     //
444     vseed->SetPoints(arrtr);
445     vseed->SetEPoints(arre);
446     //    vseed->fClusterPoints = arrcl;
447     fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
448     delete arrtr;
449     delete arre;    
450     fTreeDebug = new TTree("trackDebug","trackDebug");
451     TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
452     fTreeDebug->Branch("debug",&arrd,32000,99);
453   }
454
455
456   //set ESD event  
457   fEvent  = event;  
458 }
459
460 void AliTPCtrackerMI::FillESD(TObjArray* arr)
461 {
462   //
463   //
464   //fill esds using updated tracks
465   if (fEvent){
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);    
472       if (!pt) continue; 
473       pt->UpdatePoints();
474       //      pt->PropagateTo(fParam->GetInnerRadiusLow());
475       if (pt->GetKinkIndex(0)<=0){  //don't propagate daughter tracks 
476         pt->PropagateTo(fParam->GetInnerRadiusLow());
477       }
478  
479       if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
480         AliESDtrack iotrack;
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);
488         continue;
489       }
490        
491       if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
492         AliESDtrack iotrack;
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);
500         continue;
501       } 
502       //
503       // short tracks  - maybe decays
504
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)){
509           AliESDtrack iotrack;
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);
517           continue;
518         }
519       }       
520       
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;
526         //
527         AliESDtrack iotrack;
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);
535         continue;
536       }   
537       // short tracks  - secondaties
538       //
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){
543           AliESDtrack iotrack;
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);
551           continue;
552         }
553       }       
554       
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;
562         //
563         AliESDtrack iotrack;
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);
571         continue;
572       }   
573     }
574   }
575   printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
576 }
577
578 void AliTPCtrackerMI::WriteTracks(TTree * tree)
579 {
580   //
581   // write tracks from seed array to selected tree
582   //
583   fOutput  = tree;
584   if (fOutput){
585     AliTPCtrack *iotrack= new AliTPCtrack;
586     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
587   }
588   WriteTracks();
589 }
590
591 void AliTPCtrackerMI::WriteTracks()
592 {
593   //
594   // write tracks to the given output tree -
595   // output specified with SetIO routine
596   if (!fSeeds)  return;
597   if (!fOutput){
598     SetIO();
599   }
600
601   if (fOutput){
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;      
607     //}    
608     //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
609     TBranch * br = fOutput->GetBranch("tracks");
610     br->SetAddress(&iotrack);
611     //
612     for (Int_t i=0; i<nseed; i++) {
613       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);    
614       if (!pt) continue;    
615       AliTPCtrack * track = new AliTPCtrack(*pt);
616       iotrack = track;
617       pt->SetLab2(i); 
618       //      br->SetAddress(&iotrack);
619       fOutput->Fill();
620       delete track;
621       iotrack =0;
622     }
623     //fOutput->GetDirectory()->cd();
624     //fOutput->Write();
625   }
626   // delete iotrack;
627   //
628   if (fSeedTree){
629     //write the full seed information if specified in debug mode
630       
631     AliTPCseed * vseed = new AliTPCseed;
632     //
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);
638     //
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");
644     
645     Int_t nseed=fSeeds->GetEntriesFast();
646     
647     for (Int_t i=0; i<nseed; i++) {
648       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);    
649       if (!pt) continue;     
650       pt->SetPoints(arrtr);
651       //      pt->fClusterPoints = arrcl;
652       pt->SetEPoints(arre);
653       pt->RebuildSeed();
654       vseed = pt;
655       brseed->SetAddress(&vseed);
656       fSeedTree->Fill();
657       pt->SetPoints(0);
658       pt->SetEPoints(0);
659       //      pt->fClusterPoints = 0;
660     }
661     fSeedTree->Write();
662     if (fTreeDebug) fTreeDebug->Write();
663   }
664
665 }
666   
667
668
669
670 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
671   //
672   //
673   //seed->SetErrorY2(0.1);
674   //return 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];
684   //
685   static Float_t gcor01[500];
686   static Float_t gcor02[500];
687   static Float_t gcorp[500];
688   //
689
690   //
691   if (ginit==kFALSE){
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);
697     }
698
699     //
700     for (Int_t i=3;i<10000;i++){
701       //
702       //
703       // inner sector
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.;      
713       //
714       // outer short
715       padlength =1.;
716       gnoise2   = 0.0004/padlength;
717       nel       = 0.3*amp;
718       nprim     = 0.133*amp;
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.;
723       //
724       //
725       // outer long
726       padlength =1.5;
727       gnoise3   = 0.0004/padlength;
728       nel       = 0.3*amp;
729       nprim     = 0.133*amp;
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.;
734       //
735     }
736     ginit = kTRUE;
737   }
738   //
739   //
740   //
741   Int_t amp = int(TMath::Abs(cl->GetQ()));  
742   if (amp>9999) {
743     seed->SetErrorY2(1.);
744     return 1.;
745   }
746   Float_t snoise2;
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); 
752   //
753   //cluster "quality"
754   Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
755   Float_t res;
756   //
757   if (fSectors==fInnerSec){
758     snoise2 = gnoise1;
759     res     = ggg1[amp]*z+glandau1[amp]*angle2;     
760     if (ctype==0) res *= gcor01[rsigmay];
761     if ((ctype>0)){
762       res+=0.002;
763       res*= gcorp[rsigmay];
764     }
765   }
766   else {
767     if (padlength<1.1){
768       snoise2 = gnoise2;
769       res     = ggg2[amp]*z+glandau2[amp]*angle2; 
770       if (ctype==0) res *= gcor02[rsigmay];      
771       if ((ctype>0)){
772         res+=0.002;
773         res*= gcorp[rsigmay];
774       }
775     }
776     else{
777       snoise2 = gnoise3;      
778       res     = ggg3[amp]*z+glandau3[amp]*angle2; 
779       if (ctype==0) res *= gcor02[rsigmay];
780       if ((ctype>0)){
781         res+=0.002;
782         res*= gcorp[rsigmay];
783       }
784     }
785   }  
786
787   if (ctype<0){
788     res+=0.005;
789     res*=2.4;  // overestimate error 2 times
790   }
791   res+= snoise2;
792  
793   if (res<2*snoise2)
794     res = 2*snoise2;
795   
796   seed->SetErrorY2(res);
797   return res;
798
799
800 }
801
802
803
804 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
805   //
806   //
807   //seed->SetErrorY2(0.1);
808   //return 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];
818   //
819   static Float_t gcor01[1000];
820   static Float_t gcor02[1000];
821   static Float_t gcorp[1000];
822   //
823
824   //
825   if (ginit==kFALSE){
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);
831     }
832
833     //
834     for (Int_t i=3;i<10000;i++){
835       //
836       //
837       // inner sector
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.;      
847       //
848       // outer short
849       padlength =1.;
850       gnoise2   = 0.0004/padlength;
851       nel       = 0.3*amp;
852       nprim     = 0.133*amp;
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.;
857       //
858       //
859       // outer long
860       padlength =1.5;
861       gnoise3   = 0.0004/padlength;
862       nel       = 0.3*amp;
863       nprim     = 0.133*amp;
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.;
868       //
869     }
870     ginit = kTRUE;
871   }
872   //
873   //
874   //
875   Int_t amp = int(TMath::Abs(cl->GetQ()));  
876   if (amp>9999) {
877     seed->SetErrorY2(1.);
878     return 1.;
879   }
880   Float_t snoise2;
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());
884   //
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)); 
888   //
889   //cluster "quality"
890   Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
891   Float_t res;
892   //
893   if (fSectors==fInnerSec){
894     snoise2 = gnoise1;
895     res     = ggg1[amp]*z+glandau1[amp]*angle2;     
896     if (ctype==0) res *= gcor01[rsigmaz];
897     if ((ctype>0)){
898       res+=0.002;
899       res*= gcorp[rsigmaz];
900     }
901   }
902   else {
903     if (padlength<1.1){
904       snoise2 = gnoise2;
905       res     = ggg2[amp]*z+glandau2[amp]*angle2; 
906       if (ctype==0) res *= gcor02[rsigmaz];      
907       if ((ctype>0)){
908         res+=0.002;
909         res*= gcorp[rsigmaz];
910       }
911     }
912     else{
913       snoise2 = gnoise3;      
914       res     = ggg3[amp]*z+glandau3[amp]*angle2; 
915       if (ctype==0) res *= gcor02[rsigmaz];
916       if ((ctype>0)){
917         res+=0.002;
918         res*= gcorp[rsigmaz];
919       }
920     }
921   }  
922
923   if (ctype<0){
924     res+=0.002;
925     res*=1.3;
926   }
927   if ((ctype<0) &&amp<70){
928     res+=0.002;
929     res*=1.3;  
930   }
931   res += snoise2;
932   if (res<2*snoise2)
933      res = 2*snoise2;
934   if (res>3) res =3;
935   seed->SetErrorZ2(res);
936   return res;
937 }
938
939
940
941 /*
942 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
943   //
944   //
945   //seed->SetErrorZ2(0.1);
946   //return 0.1;
947
948   Float_t snoise2;
949   Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
950   //
951   Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
952   Int_t ctype = cl->GetType();
953   Float_t amp = TMath::Abs(cl->GetQ());
954   
955   Float_t nel;
956   Float_t nprim;
957   //
958   Float_t landau=2 ;    //landau fluctuation part
959   Float_t gg=2;         // gg fluctuation part
960   Float_t padlength= GetPadPitchLength(seed->GetX());
961  
962   if (fSectors==fInnerSec){
963     snoise2 = 0.0004/padlength;
964     nel     = 0.268*amp;
965     nprim   = 0.155*amp;
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;
969   }
970   else {
971     snoise2 = 0.0004/padlength;
972     nel     = 0.3*amp;
973     nprim   = 0.133*amp;
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;
977   }
978   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
979
980   //
981   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
982   angle2 = TMath::Sqrt((1-angle2));
983   if (angle2<0.6) angle2 = 0.6;
984   //angle2 = 1;
985
986   Float_t angle = seed->GetTgl()/angle2;
987   Float_t angular = landau*angle*angle*padlength*padlength/12.;
988   Float_t res = sdiff + angular;
989
990   
991   if ((ctype==0) && (fSectors ==fOuterSec))
992     res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
993
994   if ((ctype==0) && (fSectors ==fInnerSec))
995     res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
996   
997   if ((ctype>0)){
998     res+=0.005;
999     res*= TMath::Power(rsigmaz+0.5,1.5);  //0.31+0.147*ctype;
1000   }
1001   if (ctype<0){
1002     res+=0.002;
1003     res*=1.3;
1004   }
1005   if ((ctype<0) &&amp<70){
1006     res+=0.002;
1007     res*=1.3;  
1008   }
1009   res += snoise2;
1010   if (res<2*snoise2)
1011      res = 2*snoise2;
1012
1013   seed->SetErrorZ2(res);
1014   return res;
1015 }
1016 */
1017
1018
1019
1020
1021 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
1022 {
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());
1027   
1028   if (y > ymax) {
1029     seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
1030     if (!seed->Rotate(fSectors->GetAlpha())) 
1031       return;
1032   } else if (y <-ymax) {
1033     seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
1034     if (!seed->Rotate(-fSectors->GetAlpha())) 
1035       return;
1036   }   
1037
1038 }
1039
1040
1041
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) 
1046 {
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));
1055
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); 
1059 }
1060
1061
1062
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) 
1067 {
1068   //-----------------------------------------------------------------
1069   // Initial approximation of the track curvature
1070   //-----------------------------------------------------------------
1071   x3 -=x1;
1072   x2 -=x1;
1073   y3 -=y1;
1074   y2 -=y1;
1075   //  
1076   Double_t det = x3*y2-x2*y3;
1077   if (det==0) {
1078     return 100;
1079   }
1080   //
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);
1085   if (det<0) c2*=-1;
1086   return c2;
1087 }
1088
1089
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) 
1093 {
1094   //-----------------------------------------------------------------
1095   // Initial approximation of the track curvature
1096   //-----------------------------------------------------------------
1097   x3 -=x1;
1098   x2 -=x1;
1099   y3 -=y1;
1100   y2 -=y1;
1101   //  
1102   Double_t det = x3*y2-x2*y3;
1103   if (det==0) {
1104     return 100;
1105   }
1106   //
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);
1111   if (det<0) c2*=-1;
1112   x0+=x1;
1113   x0*=c2;  
1114   return x0;
1115 }
1116
1117
1118
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) 
1123 {
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));
1132
1133   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1134   
1135   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1136 }
1137
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) 
1142 {
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));
1147 }
1148
1149
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) 
1153 {
1154   //-----------------------------------------------------------------
1155   // Initial approximation of the tangent of the track dip angle
1156   //-----------------------------------------------------------------
1157
1158   //  Double_t angle1;
1159   
1160   //angle1    =  (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1161   //
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);
1167
1168   angle2  = (z1-z2)*c/(angle2*2.);
1169   return angle2;
1170 }
1171
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   //-----------------------------------------------------------------
1176   
1177   Double_t dx=x2-x1;
1178
1179   if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {   
1180     return kFALSE;
1181   }
1182
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);  
1185   y = x[0];
1186   z = x[1];
1187   
1188   Double_t dy = dx*(c1+c2)/(r1+r2);
1189   Double_t dz = 0;
1190   //
1191   Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1192   
1193   if (TMath::Abs(delta)>0.01){
1194     dz = x[3]*TMath::ASin(delta)/x[4];
1195   }else{
1196     dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1197   }
1198   
1199   //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1200
1201   y+=dy;
1202   z+=dz;
1203   
1204   return kTRUE;  
1205 }
1206
1207 Int_t  AliTPCtrackerMI::LoadClusters (TTree *tree)
1208 {
1209   //
1210   //
1211   fInput = tree;
1212   return LoadClusters();
1213 }
1214
1215 Int_t  AliTPCtrackerMI::LoadClusters()
1216 {
1217   //
1218   // load clusters to the memory
1219   AliTPCClustersRow *clrow= new AliTPCClustersRow;
1220   clrow->SetClass("AliTPCclusterMI");
1221   clrow->SetArray(0);
1222   clrow->GetArray()->ExpandCreateFast(10000);
1223   //
1224   //  TTree * tree = fClustersArray.GetTree();
1225
1226   TTree * tree = fInput;
1227   TBranch * br = tree->GetBranch("Segment");
1228   br->SetAddress(&clrow);
1229   //
1230   Int_t j=Int_t(tree->GetEntries());
1231   for (Int_t i=0; i<j; i++) {
1232     br->GetEntry(i);
1233     //  
1234     Int_t sec,row;
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)));
1238     }
1239     //
1240     AliTPCRow * tpcrow=0;
1241     Int_t left=0;
1242     if (sec<fkNIS*2){
1243       tpcrow = &(fInnerSec[sec%fkNIS][row]);    
1244       left = sec/fkNIS;
1245     }
1246     else{
1247       tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1248       left = (sec-fkNIS*2)/fkNOS;
1249     }
1250     if (left ==0){
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)));
1255     }
1256     if (left ==1){
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)));
1261     }
1262   }
1263   //
1264   delete clrow;
1265   LoadOuterSectors();
1266   LoadInnerSectors();
1267   return 0;
1268 }
1269
1270
1271 void AliTPCtrackerMI::UnloadClusters()
1272 {
1273   //
1274   // unload clusters from the memory
1275   //
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]);
1280       //      if (tpcrow){
1281       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
1282       //        if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
1283       //}
1284       tpcrow->ResetClusters();
1285     }
1286   //
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]);
1291       //if (tpcrow){
1292       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
1293       //if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
1294       //}
1295       tpcrow->ResetClusters();
1296     }
1297
1298   return ;
1299 }
1300
1301 void   AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1302   //
1303   //
1304   //
1305   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1306   if (!transform) {
1307     AliFatal("Tranformations not in calibDB");
1308   }
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){
1314       x[1]*=-1;
1315     }
1316   }
1317
1318   //
1319   // in debug mode  check the transformation
1320   //
1321   if (AliTPCReconstructor::StreamLevel()>1) {
1322     TTreeSRedirector &cstream = *fDebugStreamer;
1323     cstream<<"Transform"<<
1324       "x0="<<x[0]<<
1325       "x1="<<x[1]<<
1326       "x2="<<x[2]<<
1327       "Cl.="<<cluster<<
1328       "\n"; 
1329   }
1330   cluster->SetX(x[0]);
1331   cluster->SetY(x[1]);
1332   cluster->SetZ(x[2]);
1333   // The old stuff:
1334
1335   //
1336   // 
1337   //
1338   //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1339   TGeoHMatrix  *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1340   //TGeoHMatrix  mat;
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);
1344   else{
1345     // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1346   }
1347   cluster->SetX(posC[0]);
1348   cluster->SetY(posC[1]);
1349   cluster->SetZ(posC[2]);
1350 }
1351
1352 //_____________________________________________________________________________
1353 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1354   //-----------------------------------------------------------------
1355   // This function fills outer TPC sectors with clusters.
1356   //-----------------------------------------------------------------
1357   Int_t nrows = fOuterSec->GetNRows();
1358   UInt_t index=0;
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;
1363       //left
1364       Int_t ncl = tpcrow->GetN1();
1365       while (ncl--) {
1366         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1367         index=(((sec2<<8)+row)<<16)+ncl;
1368         tpcrow->InsertCluster(c,index);
1369       }
1370       //right
1371       ncl = tpcrow->GetN2();
1372       while (ncl--) {
1373         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1374         index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1375         tpcrow->InsertCluster(c,index);
1376       }
1377       //
1378       // write indexes for fast acces
1379       //
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
1385       }
1386       Int_t last = 0;
1387       for (Int_t i=0;i<510;i++){
1388         if (tpcrow->GetFastCluster(i)<0)
1389           tpcrow->SetFastCluster(i,last);
1390         else
1391           last = tpcrow->GetFastCluster(i);
1392       }
1393     }  
1394   fN=fkNOS;
1395   fSectors=fOuterSec;
1396   return 0;
1397 }
1398
1399
1400 //_____________________________________________________________________________
1401 Int_t  AliTPCtrackerMI::LoadInnerSectors() {
1402   //-----------------------------------------------------------------
1403   // This function fills inner TPC sectors with clusters.
1404   //-----------------------------------------------------------------
1405   Int_t nrows = fInnerSec->GetNRows();
1406   UInt_t index=0;
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]);
1410       //
1411       //left
1412       Int_t ncl = tpcrow->GetN1();
1413       while (ncl--) {
1414         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1415         index=(((sec<<8)+row)<<16)+ncl;
1416         tpcrow->InsertCluster(c,index);
1417       }
1418       //right
1419       ncl = tpcrow->GetN2();
1420       while (ncl--) {
1421         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1422         index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1423         tpcrow->InsertCluster(c,index);
1424       }
1425       //
1426       // write indexes for fast acces
1427       //
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
1433       }
1434       Int_t last = 0;
1435       for (Int_t i=0;i<510;i++){
1436         if (tpcrow->GetFastCluster(i)<0)
1437           tpcrow->SetFastCluster(i,last);
1438         else
1439           last = tpcrow->GetFastCluster(i);
1440       }
1441
1442     }  
1443    
1444   fN=fkNIS;
1445   fSectors=fInnerSec;
1446   return 0;
1447 }
1448
1449
1450
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;
1460
1461   const AliTPCRow * tpcrow=0;
1462   AliTPCclusterMI * clrow =0;
1463
1464   if (sec<0 || sec>=fkNIS*4) {
1465     AliWarning(Form("Wrong sector %d",sec));
1466     return 0x0;
1467   }
1468
1469   if (sec<fkNIS*2){
1470     tpcrow = &(fInnerSec[sec%fkNIS][row]);
1471     if (tpcrow==0) return 0;
1472
1473     if (sec<fkNIS) {
1474       if (tpcrow->GetN1()<=ncl) return 0;
1475       clrow = tpcrow->GetClusters1();
1476     }
1477     else {
1478       if (tpcrow->GetN2()<=ncl) return 0;
1479       clrow = tpcrow->GetClusters2();
1480     }
1481   }
1482   else {
1483     tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1484     if (tpcrow==0) return 0;
1485
1486     if (sec-2*fkNIS<fkNOS) {
1487       if (tpcrow->GetN1()<=ncl) return 0;
1488       clrow = tpcrow->GetClusters1();
1489     }
1490     else {
1491       if (tpcrow->GetN2()<=ncl) return 0;
1492       clrow = tpcrow->GetClusters2();
1493     }
1494   }
1495
1496   return &(clrow[ncl]);      
1497   
1498 }
1499
1500
1501
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   //-----------------------------------------------------------------
1506   //
1507   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1508   AliTPCclusterMI *cl=0;
1509   Int_t tpcindex= t.GetClusterIndex2(nr);
1510   //
1511   // update current shape info every 5 pad-row
1512   //  if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1513     GetShape(&t,nr);    
1514     //}
1515   //  
1516   if (fIteration>0 && tpcindex>=-1){  //if we have already clusters 
1517     //        
1518     if (tpcindex==-1) return 0; //track in dead zone
1519     if (tpcindex>0){     //
1520       cl = t.GetClusterPointer(nr);
1521       if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1522       t.SetCurrentClusterIndex1(tpcindex); 
1523     }
1524     if (cl){      
1525       Int_t relativesector = ((tpcindex&0xff000000)>>24)%18;  // if previously accepted cluster in different sector
1526       Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1527       //
1528       if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1529       if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1530       
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;      
1535       }
1536       if (!t.PropagateTo(x)) return 0;
1537       //
1538       t.SetCurrentCluster(cl); 
1539       t.SetRow(nr);
1540       Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1541       if ((tpcindex&0x8000)==0) accept =0;
1542       if (accept<3) { 
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); 
1549         }
1550         t.SetNFoundable(t.GetNFoundable()+1);
1551         UpdateTrack(&t,accept);
1552         return 1;
1553       }    
1554     }
1555   }
1556   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;  // cut on angle
1557   if (fIteration>1){
1558     // not look for new cluster during refitting    
1559     t.SetNFoundable(t.GetNFoundable()+1);
1560     return 0;
1561   }
1562   //
1563   UInt_t index=0;
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){
1567     if (y > ymax) {
1568       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1569       if (!t.Rotate(fSectors->GetAlpha())) 
1570         return 0;
1571     } else if (y <-ymax) {
1572       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1573       if (!t.Rotate(-fSectors->GetAlpha())) 
1574         return 0;
1575     }
1576     //return 1;
1577   }
1578   //
1579   if (!t.PropagateTo(x)) {
1580     if (fIteration==0) t.SetRemoval(10);
1581     return 0;
1582   }
1583   y=t.GetY(); 
1584   Double_t z=t.GetZ();
1585   //
1586
1587   if (!IsActive(t.GetRelativeSector(),nr)) {
1588     t.SetInDead(kTRUE);
1589     t.SetClusterIndex2(nr,-1); 
1590     return 0;
1591   }
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;
1595   
1596   if (!isActive || !isActive2) return 0;
1597
1598   const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1599   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1600   Double_t  roady  =1.;
1601   Double_t  roadz = 1.;
1602   //
1603   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1604     t.SetInDead(kTRUE);
1605     t.SetClusterIndex2(nr,-1); 
1606     return 0;
1607   } 
1608   else
1609     {
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);
1612       else
1613         return 0;
1614     }   
1615   //calculate 
1616   if (krow) {
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));       
1620   }  
1621   if (cl) {
1622     t.SetCurrentCluster(cl); 
1623     t.SetRow(nr);
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); 
1631     }
1632     /*    
1633     if (t.fCurrentCluster->IsUsed(10)){
1634       //
1635       //     
1636
1637       t.fNShared++;
1638       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1639         t.fRemoval =10;
1640         return 0;
1641       }
1642     }
1643     */
1644     if (accept<3) UpdateTrack(&t,accept);
1645
1646   } else {  
1647     if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1648     
1649   }
1650   return 1;
1651 }
1652
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   //-----------------------------------------------------------------
1657   //
1658   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1659   Double_t y,z; 
1660   if (!t.GetProlongation(x,y,z)) {
1661     t.SetRemoval(10);
1662     return 0;
1663   }
1664   //
1665   //
1666   if (TMath::Abs(y)>ymax){
1667     
1668     if (y > ymax) {
1669       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1670       if (!t.Rotate(fSectors->GetAlpha())) 
1671         return 0;
1672     } else if (y <-ymax) {
1673       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1674       if (!t.Rotate(-fSectors->GetAlpha())) 
1675         return 0;
1676     }
1677     if (!t.PropagateTo(x)) {
1678       return 0;
1679     } 
1680     t.GetProlongation(x,y,z);
1681   }
1682   //
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);
1687     GetShape(&t,nr);
1688   }
1689   //  
1690   AliTPCclusterMI *cl=0;
1691   UInt_t index=0;
1692   
1693   
1694   //Int_t nr2 = nr;
1695   const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1696   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1697   Double_t  roady  =1.;
1698   Double_t  roadz = 1.;
1699   //
1700   Int_t row = nr;
1701   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1702     t.SetInDead(kTRUE);
1703     t.SetClusterIndex2(row,-1); 
1704     return 0;
1705   } 
1706   else
1707     {
1708       if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1709     }   
1710   //calculate 
1711   
1712   if ((cl==0)&&(krow)) {
1713     //    cl = krow.FindNearest2(y+10,z,roady,roadz,index);    
1714     cl = krow.FindNearest2(y,z,roady,roadz,index);    
1715
1716     if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));       
1717   }  
1718
1719   if (cl) {
1720     t.SetCurrentCluster(cl); 
1721     //    Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);        
1722     //if (accept<3){
1723       t.SetClusterIndex2(row,index);
1724       t.SetClusterPointer(row, cl);
1725       //}
1726   }
1727   return 1;
1728 }
1729
1730
1731
1732 //_________________________________________________________________________
1733 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1734 {
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;
1741   Float_t xyz[3];
1742   //  xyz[0] = fParam->GetPadRowRadii(sector,row);
1743   xyz[0] = cl->GetX();
1744   xyz[1] = cl->GetY();
1745   xyz[2] = cl->GetZ();
1746   Float_t sin,cos;
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];
1750   Float_t cov[6];
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;
1757   cov[2] = 0.;
1758   cov[3] = cos*cos*sigmaY2;
1759   cov[4] = 0.;
1760   cov[5] = sigmaZ2;
1761   p.SetXYZ(x,y,xyz[2],cov);
1762   AliGeomManager::ELayerID iLayer;
1763   Int_t idet;
1764   if (sector < fParam->GetNInnerSector()) {
1765     iLayer = AliGeomManager::kTPC1;
1766     idet = sector;
1767   }
1768   else {
1769     iLayer = AliGeomManager::kTPC2;
1770     idet = sector - fParam->GetNInnerSector();
1771   }
1772   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1773   p.SetVolumeID(volid);
1774   return kTRUE;
1775 }
1776
1777
1778
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);   
1785    
1786   Double_t xt=t.GetX();
1787   Int_t     row = GetRowNumber(xt)-1; 
1788   Double_t  ymax= GetMaxY(nr);
1789
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) {
1792 //     t.fRemoval =10;
1793 //     return 0;  // not prolongate strongly inclined tracks
1794 //   } 
1795 //   if (TMath::Abs(t.GetSnp())>0.95) {
1796 //     t.fRemoval =10;
1797 //     return 0;  // not prolongate strongly inclined tracks
1798 //   }// patch 28 fev 06
1799
1800   Double_t x= GetXrow(nr);
1801   Double_t y,z;
1802   //t.PropagateTo(x+0.02);
1803   //t.PropagateTo(x+0.01);
1804   if (!t.PropagateTo(x)){
1805     return 0;
1806   }
1807   //
1808   y=t.GetY();
1809   z=t.GetZ();
1810
1811   if (TMath::Abs(y)>ymax){
1812     if (y > ymax) {
1813       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1814       if (!t.Rotate(fSectors->GetAlpha())) 
1815         return 0;
1816     } else if (y <-ymax) {
1817       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1818       if (!t.Rotate(-fSectors->GetAlpha())) 
1819         return 0;
1820     }
1821     //    if (!t.PropagateTo(x)){
1822     //  return 0;
1823     //}
1824     return 1;
1825     //y = t.GetY();    
1826   }
1827   //
1828   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1829
1830   if (!IsActive(t.GetRelativeSector(),nr)) {
1831     t.SetInDead(kTRUE);
1832     t.SetClusterIndex2(nr,-1); 
1833     return 0;
1834   }
1835   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1836
1837   AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1838
1839   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1840     t.SetInDead(kTRUE);
1841     t.SetClusterIndex2(nr,-1); 
1842     return 0;
1843   } 
1844   else
1845     {
1846       if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1847       else
1848         return 0;      
1849     }
1850
1851   // update current
1852   if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1853     //    t.fCurrentSigmaY = GetSigmaY(&t);
1854     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1855     GetShape(&t,nr);
1856   }
1857     
1858   AliTPCclusterMI *cl=0;
1859   Int_t index=0;
1860   //
1861   Double_t roady = 1.;
1862   Double_t roadz = 1.;
1863   //
1864
1865   if (!cl){
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);
1871       if (cl) {
1872         t.SetCurrentCluster(cl);
1873         return 1;
1874       }
1875     }
1876   }
1877
1878   //  if (index<0) return 0;
1879   UInt_t uindex = TMath::Abs(index);
1880
1881   if (krow) {    
1882     //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);      
1883     cl = krow.FindNearest2(y,z,roady,roadz,uindex);      
1884   }
1885
1886   if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));   
1887   t.SetCurrentCluster(cl);
1888
1889   return 1;
1890 }
1891
1892
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   //-----------------------------------------------------------------
1897
1898   //update error according neighborhoud
1899
1900   if (t.GetCurrentCluster()) {
1901     t.SetRow(nr); 
1902     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1903     
1904     if (t.GetCurrentCluster()->IsUsed(10)){
1905       //
1906       //
1907       //  t.fErrorZ2*=2;
1908       //  t.fErrorY2*=2;
1909       t.SetNShared(t.GetNShared()+1);
1910       if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1911         t.SetRemoval(10);
1912         return 0;
1913       }
1914     }   
1915     if (fIteration>0) accept = 0;
1916     if (accept<3)  UpdateTrack(&t,accept);  
1917  
1918   } else {
1919     if (fIteration==0){
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);      
1922
1923       if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1924     }
1925   }
1926   return 1;
1927 }
1928
1929
1930
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();
1937   //
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();  
1941   //
1942   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1943     
1944   Int_t first = GetRowNumber(xt)-1;
1945   for (Int_t nr= first; nr>=rf; nr-=step) {
1946     // update kink info
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;
1952         //
1953         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1954         if (!kink){
1955           printf("PROBLEM\n");
1956         }
1957         else{
1958           Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1959           if (kinkrow==nr){
1960             AliExternalTrackParam paramd(t);
1961             kink->SetDaughter(paramd);
1962             kink->SetStatus(2,5);
1963             kink->Update();
1964           }
1965         }
1966       }
1967     }
1968
1969     if (nr==80) t.UpdateReference();
1970     if (nr<fInnerSec->GetNRows()) 
1971       fSectors = fInnerSec;
1972     else
1973       fSectors = fOuterSec;
1974     if (FollowToNext(t,nr)==0) 
1975       if (!t.IsActive()) 
1976         return 0;
1977     
1978   }   
1979   return 1;
1980 }
1981
1982
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();
1989   //
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);
1994     
1995   for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1996     
1997     if (FollowToNextFast(t,nr)==0) 
1998       if (!t.IsActive()) return 0;
1999     
2000   }   
2001   return 1;
2002 }
2003
2004
2005
2006
2007
2008 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
2009   //-----------------------------------------------------------------
2010   // This function tries to find a track prolongation.
2011   //-----------------------------------------------------------------
2012   //
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);
2018     
2019   Int_t first = t.GetFirstPoint();
2020   if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
2021   //
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);
2032         if (!kink){
2033           printf("PROBLEM\n");
2034         }
2035         else{
2036           Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2037           if (kinkrow==nr){
2038             AliExternalTrackParam paramm(t);
2039             kink->SetMother(paramm);
2040             kink->SetStatus(2,1);
2041             kink->Update();
2042           }
2043         }
2044       }      
2045     }
2046     //
2047     if (nr<fInnerSec->GetNRows()) 
2048       fSectors = fInnerSec;
2049     else
2050       fSectors = fOuterSec;
2051
2052     FollowToNext(t,nr);                                                             
2053   }   
2054   return 1;
2055 }
2056
2057
2058
2059
2060    
2061 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2062 {
2063   //
2064   //
2065   sum1=0;
2066   sum2=0;
2067   Int_t sum=0;
2068   //
2069   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2070   dz2*=dz2;  
2071
2072   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2073   dy2*=dy2;
2074   Float_t distance = TMath::Sqrt(dz2+dy2);
2075   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
2076  
2077   //  Int_t offset =0;
2078   Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2079   Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2080   if (lastpoint>160) 
2081     lastpoint =160;
2082   if (firstpoint<0) 
2083     firstpoint = 0;
2084   if (firstpoint>lastpoint) {
2085     firstpoint =lastpoint;
2086     //    lastpoint  =160;
2087   }
2088     
2089   
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) {
2094       sum++;
2095     }
2096   }
2097   if (sum<5) return 0;
2098
2099   Float_t summin = TMath::Min(sum1+1,sum2+1);
2100   Float_t ratio = (sum+1)/Float_t(summin);
2101   return ratio;
2102 }
2103
2104 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2105 {
2106   //
2107   //
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));
2112   
2113   //
2114   Int_t sumshared=0;
2115   //
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;
2120   //
2121   //  if (firstpoint>=lastpoint-5) return;;
2122
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) {
2126       sumshared++;
2127       s1->SetSharedMapBit(i, kTRUE);
2128       s2->SetSharedMapBit(i, kTRUE);
2129     }
2130     if (s1->GetClusterIndex2(i)>0)
2131       s1->SetClusterMapBit(i, kTRUE);
2132   }
2133   if (sumshared>cutN0){
2134     // sign clusters
2135     //
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);
2144         }       
2145       }
2146     }
2147   }
2148   //  
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());
2155         break;
2156       } 
2157     }
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());
2163         break;
2164       } 
2165     }    
2166   }
2167 }
2168
2169 void  AliTPCtrackerMI::SignShared(TObjArray * arr)
2170 {
2171   //
2172   //sort trackss according sectors
2173   //  
2174   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2175     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2176     if (!pt) continue;
2177     //if (pt) RotateToLocal(pt);
2178     pt->SetSort(0);
2179   }
2180   arr->UnSort();
2181   arr->Sort();  // sorting according relative sectors
2182   arr->Expand(arr->GetEntries());
2183   //
2184   //
2185   Int_t nseed=arr->GetEntriesFast();
2186   for (Int_t i=0; i<nseed; i++) {
2187     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2188     if (!pt) continue;
2189     for (Int_t j=0;j<=12;j++){
2190       pt->SetOverlapLabel(j,0);
2191     }
2192   }
2193   for (Int_t i=0; i<nseed; i++) {
2194     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2195     if (!pt) continue;
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;
2200       //      if (pt2){
2201       if (pt2->GetRemoval()<=10) {
2202         //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2203         SignShared(pt,pt2);
2204       }
2205     }  
2206   }
2207 }
2208
2209
2210 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2211 {
2212   //
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);    
2217     if (!pt) {
2218       continue;
2219     }
2220     pt->SetSort(mode);
2221   }
2222   arr->UnSort();
2223   arr->Sort();
2224 }
2225
2226
2227 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t minimal)
2228 {
2229   //
2230   // Loop over all tracks and remove overlaped tracks (with lower quality)
2231   // Algorithm:
2232   //    1. Unsign clusters
2233   //    2. Sort tracks according quality
2234   //       Quality is defined by the number of cluster between first and last points 
2235   //       
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
2240   //
2241   //Called in - AliTPCtrackerMI::Clusters2Tracks()
2242   //          - AliTPCtrackerMI::PropagateBack()
2243   //          - AliTPCtrackerMI::RefitInward()
2244   //
2245
2246  
2247   UnsignClusters();
2248   //
2249   Int_t nseed = arr->GetEntriesFast();  
2250   Float_t * quality = new Float_t[nseed];
2251   Int_t   * indexes = new Int_t[nseed];
2252   Int_t good =0;
2253   //
2254   //
2255   for (Int_t i=0; i<nseed; i++) {
2256     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2257     if (!pt){
2258       quality[i]=-1;
2259       continue;
2260     }
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); 
2267   }
2268   TMath::Sort(nseed,quality,indexes);
2269   //
2270   //
2271   for (Int_t itrack=0; itrack<nseed; itrack++) {
2272     Int_t trackindex = indexes[itrack];
2273     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);   
2274     if (!pt) continue;
2275     //
2276     if (quality[trackindex]<0){
2277       if (pt) {
2278         delete arr->RemoveAt(trackindex);
2279       }
2280       else{
2281         arr->RemoveAt(trackindex);
2282       }
2283       continue;
2284     }
2285     //
2286     //
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;
2290     //
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;
2295     if (pt->GetESD()){
2296       Int_t dummy[12];
2297       if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2298     }
2299     if (!itsgold){
2300       //
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);
2304         continue;
2305       }      
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);
2309         continue;
2310       }
2311     }
2312
2313     good++;
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 ? 
2318     //
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);        
2324       if (!c) continue;
2325       c->Use(10);  
2326     }    
2327   }
2328   fNtracks = good;
2329   if (fDebug>0){
2330     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2331   }
2332   delete []quality;
2333   delete []indexes;
2334 }
2335
2336 void AliTPCtrackerMI::UnsignClusters() 
2337 {
2338   //
2339   // loop over all clusters and unsign them
2340   //
2341   
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))         
2347         cl[icl].Use(-1);
2348       cl = fInnerSec[sec][row].GetClusters2();
2349       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2350         //if (cl[icl].IsUsed(10))       
2351           cl[icl].Use(-1);      
2352     }
2353   }
2354   
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))       
2360           cl[icl].Use(-1);
2361       cl = fOuterSec[sec][row].GetClusters2();
2362       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2363         //if (cl[icl].IsUsed(10))       
2364         cl[icl].Use(-1);      
2365     }
2366   }
2367   
2368 }
2369
2370
2371
2372 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2373 {
2374   //
2375   //sign clusters to be "used"
2376   //
2377   // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2378   // loop over "primaries"
2379   
2380   Float_t sumdens=0;
2381   Float_t sumdens2=0;
2382   Float_t sumn   =0;
2383   Float_t sumn2  =0;
2384   Float_t sumchi =0;
2385   Float_t sumchi2 =0;
2386
2387   Float_t sum    =0;
2388
2389   TStopwatch timer;
2390   timer.Start();
2391
2392   Int_t nseed = arr->GetEntriesFast();
2393   for (Int_t i=0; i<nseed; i++) {
2394     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2395     if (!pt) {
2396       continue;
2397     }    
2398     if (!(pt->IsActive())) continue;
2399     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2400     if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2401       sumdens += dens;
2402       sumdens2+= dens*dens;
2403       sumn    += pt->GetNumberOfClusters();
2404       sumn2   += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2405       Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2406       if (chi2>5) chi2=5;
2407       sumchi  +=chi2;
2408       sumchi2 +=chi2*chi2;
2409       sum++;
2410     }
2411   }
2412
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;
2419   
2420
2421   if (sum>20){
2422     mdensity = sumdens/sum;
2423     meann    = sumn/sum;
2424     meanchi  = sumchi/sum;
2425     //
2426     sdensity = sumdens2/sum-mdensity*mdensity;
2427     if (sdensity >= 0)
2428        sdensity = TMath::Sqrt(sdensity);
2429     else
2430        sdensity = 0.1;
2431     //
2432     smeann   = sumn2/sum-meann*meann;
2433     if (smeann >= 0)
2434       smeann   = TMath::Sqrt(smeann);
2435     else 
2436       smeann   = 10;
2437     //
2438     smeanchi = sumchi2/sum - meanchi*meanchi;
2439     if (smeanchi >= 0)
2440       smeanchi = TMath::Sqrt(smeanchi);
2441     else
2442       smeanchi = 0.4;
2443   }
2444
2445
2446   //REMOVE  SHORT DELTAS or tracks going out of sensitive volume of TPC
2447   //
2448   for (Int_t i=0; i<nseed; i++) {
2449     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2450     if (!pt) {
2451       continue;
2452     }
2453     if (pt->GetBSigned()) continue;
2454     if (pt->GetBConstrain()) continue;    
2455     //if (!(pt->IsActive())) continue;
2456     /*
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);
2462       }
2463       continue;
2464     }
2465     */
2466     Bool_t isok =kFALSE;
2467     if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2468       isok = kTRUE;
2469     if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2470       isok =kTRUE;
2471     if  (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2472       isok =kTRUE;
2473     if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2474       isok =kTRUE;
2475     
2476     if (isok)     
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);
2481         if (!c) continue;
2482         //if (!(c->IsUsed(10))) c->Use();  
2483         c->Use(10);  
2484       }
2485   }
2486   
2487   
2488   //
2489   Double_t maxchi  = meanchi+2.*smeanchi;
2490
2491   for (Int_t i=0; i<nseed; i++) {
2492     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2493     if (!pt) {
2494       continue;
2495     }    
2496     //if (!(pt->IsActive())) continue;
2497     if (pt->GetBSigned()) continue;
2498     Double_t chi     = pt->GetChi2()/pt->GetNumberOfClusters();
2499     if (chi>maxchi) continue;
2500
2501     Float_t bfactor=1;
2502     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2503    
2504     //sign only tracks with enoug big density at the beginning
2505     
2506     if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue; 
2507     
2508     
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);
2511    
2512     //    if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2513     if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2514       minn=0;
2515
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++) {
2520
2521         Int_t index=pt->GetClusterIndex2(i);
2522         if (index<0) continue;
2523         AliTPCclusterMI *c= pt->GetClusterPointer(i);
2524         if (!c) continue;
2525         //      if (!(c->IsUsed(10))) c->Use();  
2526         c->Use(10);  
2527       }
2528     }
2529   }
2530   //  gLastCheck = nseed;
2531   //  arr->Compress();
2532   if (fDebug>0){
2533     timer.Print();
2534   }
2535 }
2536
2537
2538 void  AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2539 {
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();  
2544   //
2545   for (Int_t i=0; i<nseed; i++) {
2546     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
2547     if (!pt) {
2548       continue;
2549     }
2550     if (!(pt->IsActive())) continue;
2551     StopNotActive(pt,row0,th0, th1,th2);
2552   }
2553 }
2554
2555
2556
2557 void  AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2558  Float_t th2) const
2559 {
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 
2563   Int_t sumgood1  = 0;
2564   Int_t sumgood2  = 0;
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) ;
2569     return;
2570   }
2571
2572   for (Int_t i=row0; i<maxindex; i++){
2573     Int_t index = seed->GetClusterIndex2(i);
2574     if (index!=-1) foundable++;
2575     //if (!c) continue;
2576     if (foundable<=30) sumgood1++;
2577     if (foundable<=50) {
2578       sumgood2++;
2579     }
2580     else{ 
2581       break;
2582     }        
2583   }
2584   if (foundable>=30.){ 
2585      if (sumgood1<(th1*30.)) seed->Desactivate(10);
2586   }
2587   if (foundable>=50)
2588     if (sumgood2<(th2*50.)) seed->Desactivate(10);
2589 }
2590
2591
2592 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2593 {
2594   //
2595   // back propagation of ESD tracks
2596   //
2597   //return 0;
2598   const Int_t kMaxFriendTracks=2000;
2599   fEvent = event;
2600   ReadSeeds(event,2);
2601   fIteration=2;
2602   //PrepareForProlongation(fSeeds,1);
2603   PropagateForward2(fSeeds);
2604   RemoveUsed2(fSeeds,0.4,0.4,20);
2605
2606   TObjArray arraySeed(fSeeds->GetEntries());
2607   for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2608     arraySeed.AddAt(fSeeds->At(i),i);    
2609   }
2610   SignShared(&arraySeed);
2611   //  FindCurling(fSeeds, event,2); // find multi found tracks
2612   FindSplitted(fSeeds, event,2); // find multi found tracks
2613
2614   Int_t ntracks=0;
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
2620
2621     seed->PropagateTo(fParam->GetInnerRadiusLow());
2622     seed->UpdatePoints();
2623     MakeBitmaps(seed);
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());
2629     //
2630     if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2631       TTreeSRedirector &cstream = *fDebugStreamer;
2632       cstream<<"Crefit"<<
2633         "Esd.="<<esd<<
2634         "Track.="<<seed<<
2635         "\n"; 
2636     }
2637
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);
2646       //
2647       // add seed to the esd track in Calib level
2648       //
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);
2653       }
2654       ntracks++;
2655     }
2656     else{
2657       //printf("problem\n");
2658     }
2659   }
2660   //FindKinks(fSeeds,event);
2661   Info("RefitInward","Number of refitted tracks %d",ntracks);
2662   fEvent =0;
2663   //WriteTracks();
2664   return 0;
2665 }
2666
2667
2668 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2669 {
2670   //
2671   // back propagation of ESD tracks
2672   //
2673
2674   fEvent = event;
2675   fIteration = 1;
2676   ReadSeeds(event,1);
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
2681
2682   //
2683   Int_t nseed = fSeeds->GetEntriesFast();
2684   Int_t ntracks=0;
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);
2701       ntracks++;
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"<<
2706           "Tr0.="<<seed<<
2707           "EventNrInFile="<<eventnumber<<
2708           "\n"; // patch 28 fev 06         
2709       }
2710     }
2711   }
2712   //FindKinks(fSeeds,event);
2713   Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2714   fEvent =0;
2715   //WriteTracks();
2716   
2717   return 0;
2718 }
2719
2720
2721 void AliTPCtrackerMI::DeleteSeeds()
2722 {
2723   //
2724   //delete Seeds
2725
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);
2730   }
2731   delete fSeeds;
2732
2733   fSeeds =0;
2734 }
2735
2736 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2737 {
2738   //
2739   //read seeds from the event
2740   
2741   Int_t nentr=event->GetNumberOfTracks();
2742   if (fDebug>0){
2743     Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2744   }
2745   if (fSeeds) 
2746     DeleteSeeds();
2747   if (!fSeeds){   
2748     fSeeds = new TObjArray(nentr);
2749   }
2750   UnsignClusters();
2751   //  Int_t ntrk=0;  
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);
2770       }
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);
2774       }
2775
2776     }
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) ) {
2780       fSeeds->AddAt(0,i);
2781       delete seed;
2782       continue;    
2783     }
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]);
2790       Double_t trdchi2=0;
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.);
2795     }
2796
2797     //
2798     //
2799     // rotate to the local coordinate system
2800     //   
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)) {
2811       delete seed;
2812       continue;
2813     }
2814     seed->SetESD(esd);
2815     // sign clusters
2816     if (esd->GetKinkIndex(0)<=0){
2817       for (Int_t irow=0;irow<160;irow++){
2818         Int_t index = seed->GetClusterIndex2(irow);    
2819         if (index>0){ 
2820           //
2821           AliTPCclusterMI * cl = GetClusterMI(index);
2822           seed->SetClusterPointer(irow,cl);
2823           if (cl){
2824             if ((index & 0x8000)==0){
2825               cl->Use(10);  // accepted cluster   
2826             }else{
2827               cl->Use(6);   // close cluster not accepted
2828             }   
2829           }else{
2830             Info("ReadSeeds","Not found cluster");
2831           }
2832         }
2833       }
2834     }
2835     fSeeds->AddAt(seed,i);
2836   }
2837 }
2838
2839
2840
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
2852   Int_t nin0  = 0;
2853   Int_t nin1  = 0;
2854   Int_t nin2  = 0;
2855   Int_t nin   = 0;
2856   Int_t nout1 = 0;
2857   Int_t nout2 = 0;
2858
2859   Double_t x[5], c[15];
2860   //  Int_t di = i1-i2;
2861   //
2862   AliTPCseed * seed = new AliTPCseed();
2863   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2864   Double_t cs=cos(alpha), sn=sin(alpha);
2865   //
2866   //  Double_t x1 =fOuterSec->GetX(i1);
2867   //Double_t xx2=fOuterSec->GetX(i2);
2868   
2869   Double_t x1 =GetXrow(i1);
2870   Double_t xx2=GetXrow(i2);
2871
2872   Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2873
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
2877   //
2878   Int_t ns =sec;   
2879
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;  
2883
2884   //
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.);
2890   }
2891   Double_t r2min = 1/(cuts[0]*cuts[0]);  //minimal square of radius given by cut
2892
2893   //  Int_t ddsec = 1;
2894   if (deltay>0) ddsec = 0; 
2895   // loop over clusters  
2896   for (Int_t is=0; is < kr1; is++) {
2897     //
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;
2901
2902     if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
2903
2904     // find possible directions    
2905     Float_t anglez = (z1-z3)/(x1-x3); 
2906     Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z      
2907     //
2908     //
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);            
2914     
2915     //
2916     // loop over 2 sectors
2917     Int_t dsec1=-ddsec;
2918     Int_t dsec2= ddsec;
2919     if (y1<0)  dsec2= 0;
2920     if (y1>0)  dsec1= 0;
2921     
2922     Double_t dddz1=0;  // direction of delta inclination in z axis
2923     Double_t dddz2=0;
2924     if ( (z1-z3)>0)
2925       dddz1 =1;    
2926     else
2927       dddz2 =1;
2928     //
2929     for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2930       Int_t sec2 = sec + dsec;
2931       // 
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);
2938
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; 
2942       //
2943       //      Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2944
2945       //
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;  
2951         //
2952         //calcutate parameters
2953         //      
2954         Double_t yy0 =  dyy0 +(kcl->GetY()-y3)*cs13r;
2955         // stright track
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;                
2962        
2963         nin0++; 
2964         //
2965         Double_t c0  = 1/TMath::Sqrt(r02);
2966         if (yy0>0) c0*=-1.;     
2967                
2968        
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);  
2973         //
2974         //
2975         Double_t z0  =  kcl->GetZ();  
2976         Double_t zzzz2    = z1-(z1-z3)*dfi1/dfi0;
2977         if (TMath::Abs(zzzz2-z0)>0.5) continue;       
2978         nin1++;              
2979         //      
2980         Double_t dip    = (z1-z0)*c0/dfi1;        
2981         Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2982         //
2983         y2 = kcl->GetY(); 
2984         if (dsec==0){
2985           x2 = xx2; 
2986           z2 = kcl->GetZ();       
2987         }
2988         else
2989           {
2990             // rotation 
2991             z2 = kcl->GetZ();  
2992             x2= xx2*cs-y2*sn*dsec;
2993             y2=+xx2*sn*dsec+y2*cs;
2994           }
2995         
2996         x[0] = y1;
2997         x[1] = z1;
2998         x[2] = x0;
2999         x[3] = dip;
3000         x[4] = c0;
3001         //
3002         //
3003         // do we have cluster at the middle ?
3004         Double_t ym,zm;
3005         GetProlongation(x1,xm,x,ym,zm);
3006         UInt_t dummy; 
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))) {        
3011             continue;
3012           }
3013         }
3014         else{     
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};
3020           //
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))) {      
3025               continue;
3026             }
3027           }
3028         }
3029        
3030
3031         Double_t dym = 0;
3032         Double_t dzm = 0;
3033         if (cm){
3034           dym = ym - cm->GetY();
3035           dzm = zm - cm->GetZ();
3036         }
3037         nin2++;
3038
3039
3040         //
3041         //
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;
3047
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;
3054         
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;
3059         
3060         c[0]=sy1;
3061         c[1]=0.;       c[2]=sz1;
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;
3068         
3069         //      if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3070         
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);
3074         
3075         track->SetIsSeeding(kTRUE);
3076         track->SetSeed1(i1);
3077         track->SetSeed2(i2);
3078         track->SetSeedType(3);
3079
3080        
3081         //if (dsec==0) {
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){
3086             seed->Reset();
3087             seed->~AliTPCseed();
3088             continue;
3089           }
3090           //}
3091         
3092         nin++;
3093         FollowProlongation(*track, i2,1);
3094         
3095         
3096         //Int_t rc = 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());
3101         
3102         if (track->GetNumberOfClusters()<(i1-i2)*0.5 || 
3103             track->GetNumberOfClusters() < track->GetNFoundable()*0.6 || 
3104             track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3105           seed->Reset();
3106           seed->~AliTPCseed();
3107           continue;
3108         }
3109         nout1++;
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); 
3126               seed->Reset();
3127               seed->~AliTPCseed();
3128               continue;         
3129             }
3130             else{
3131               seed->Reset();
3132               seed->~AliTPCseed();
3133               continue;
3134             
3135             }
3136           }
3137         }
3138       
3139         track->SetSeedType(0);
3140         arr->AddLast(track); 
3141         seed = new AliTPCseed;  
3142         nout2++;
3143         // don't consider other combinations
3144         if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3145           break;
3146       }
3147     }
3148   }
3149   if (fDebug>3){
3150     Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3151   }
3152   delete seed;
3153 }
3154
3155
3156 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
3157                                  Float_t deltay) {
3158   
3159
3160
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
3168
3169
3170   Int_t nin0  = 0;
3171   Int_t nin1  = 0;
3172   Int_t nin2  = 0;
3173   Int_t nin   = 0;
3174   Int_t nout1 = 0;
3175   Int_t nout2 = 0;
3176   Int_t nout3 =0;
3177   Double_t x[5], c[15];
3178   //
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);
3183   //
3184   //
3185
3186   // first 3 padrows
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;  
3190   //
3191   Double_t x1p = GetXrow(i1);
3192   const    AliTPCRow& kr1p=GetRow(sec,i1);
3193   //
3194   Double_t x1m = GetXrow(i1-2);
3195   const    AliTPCRow& kr1m=GetRow(sec,i1-2);
3196
3197   //
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;  
3202   //
3203   AliTPCRow&  kr3p  = GetRow((sec+fkNOS)%fkNOS,i1-6);
3204   Double_t    x3p   = GetXrow(i1-6);
3205   //
3206   AliTPCRow&  kr3m  = GetRow((sec+fkNOS)%fkNOS,i1-8);
3207   Double_t    x3m   = GetXrow(i1-8);
3208
3209   //
3210   //
3211   // middle padrow
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;  
3216   //
3217   //
3218   Double_t deltax  = x1-x3;
3219   Double_t dymax   = deltax*cuts[1];
3220   Double_t dzmax   = deltax*cuts[3];
3221   //
3222   // loop over clusters  
3223   for (Int_t is=0; is < kr1; is++) {
3224     //
3225     if (kr1[is]->IsUsed(10)) continue;
3226     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();    
3227     //
3228     if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge    
3229     // 
3230     Int_t  index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3231     Int_t  index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3232     //    
3233     Double_t y3,   z3;
3234     //
3235     //
3236     UInt_t index;
3237     for (Int_t js=index1; js < index2; js++) {
3238       const AliTPCclusterMI *kcl = kr3[js];
3239       if (kcl->IsUsed(10)) continue;
3240       y3 = kcl->GetY(); 
3241       // apply angular cuts
3242       if (TMath::Abs(y1-y3)>dymax) continue;
3243       x3 = x3; 
3244       z3 = kcl->GetZ(); 
3245       if (TMath::Abs(z1-z3)>dzmax) continue;
3246       //
3247       Double_t angley = (y1-y3)/(x1-x3);
3248       Double_t anglez = (z1-z3)/(x1-x3);
3249       //
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;
3252       //
3253       Double_t yyym = angley*(xm-x1)+y1;
3254       Double_t zzzm = anglez*(xm-x1)+z1;
3255
3256       const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3257       if (!kcm) continue;
3258       if (kcm->IsUsed(10)) continue;
3259       
3260       erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3261       errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3262       //
3263       //
3264       //
3265       Int_t used  =0;
3266       Int_t found =0;
3267       //
3268       // look around first
3269       const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3270                                                       anglez*(x1m-x1)+z1,
3271                                                       erry,errz,index);
3272       //
3273       if (kc1m){
3274         found++;
3275         if (kc1m->IsUsed(10)) used++;
3276       }
3277       const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3278                                                       anglez*(x1p-x1)+z1,
3279                                                       erry,errz,index);
3280       //
3281       if (kc1p){
3282         found++;
3283         if (kc1p->IsUsed(10)) used++;
3284       }
3285       if (used>1)  continue;
3286       if (found<1) continue; 
3287
3288       //
3289       // look around last
3290       const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3291                                                       anglez*(x3m-x3)+z3,
3292                                                       erry,errz,index);
3293       //
3294       if (kc3m){
3295         found++;
3296         if (kc3m->IsUsed(10)) used++;
3297       }
3298       else 
3299         continue;
3300       const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3301                                                       anglez*(x3p-x3)+z3,
3302                                                       erry,errz,index);
3303       //
3304       if (kc3p){
3305         found++;
3306         if (kc3p->IsUsed(10)) used++;
3307       }
3308       else 
3309         continue;
3310       if (used>1)  continue;
3311       if (found<3) continue;       
3312       //
3313       Double_t x2,y2,z2;
3314       x2 = xm;
3315       y2 = kcm->GetY();
3316       z2 = kcm->GetZ();
3317       //
3318                         
3319       x[0]=y1;
3320       x[1]=z1;
3321       x[4]=F1(x1,y1,x2,y2,x3,y3);
3322       //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3323       nin0++;
3324       //
3325       x[2]=F2(x1,y1,x2,y2,x3,y3);
3326       nin1++;
3327       //
3328       x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3329       //if (TMath::Abs(x[3]) > cuts[3]) continue;
3330       nin2++;
3331       //
3332       //
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;
3336       
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;
3343       
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;
3348       
3349       c[0]=sy1;
3350       c[1]=0.;       c[2]=sz1; 
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;
3357       
3358       //        if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3359       
3360       UInt_t index=kr1.GetIndex(is);
3361       seed->~AliTPCseed();
3362       AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3363       
3364       track->SetIsSeeding(kTRUE);
3365
3366       nin++;      
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){
3370         seed->Reset();
3371         seed->~AliTPCseed();
3372         continue;
3373       }
3374       nout1++;
3375       nout2++;  
3376       //Int_t rc = 1;
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());
3381       
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 ) {
3385         seed->Reset();
3386         seed->~AliTPCseed();
3387         continue;
3388       }
3389    
3390       {
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); 
3397         seed->Reset();
3398         seed->~AliTPCseed();
3399       }
3400       
3401    
3402       //arr->AddLast(track); 
3403       //seed = new AliTPCseed;  
3404       nout3++;
3405     }
3406   }
3407   
3408   if (fDebug>3){
3409     Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3410   }
3411   delete seed;
3412 }
3413
3414
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
3425   Int_t nin0=0;
3426   Int_t nin1=0;
3427   Int_t nin2=0;
3428   Int_t nin3=0;
3429   //  Int_t nin4=0;
3430   //Int_t nin5=0;
3431
3432   
3433
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];
3439   AliTPCRow * kr=0;
3440
3441   AliTPCpolyTrack polytrack;
3442   Int_t nclusters=fSectors[sec][row0];
3443   AliTPCseed * seed = new AliTPCseed;
3444
3445   Int_t sumused=0;
3446   Int_t cused=0;
3447   Int_t cnused=0;
3448   for (Int_t is=0; is < nclusters; is++) {  //LOOP over clusters
3449     Int_t nfound =0;
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];
3455       
3456       if (cl->IsUsed(10)) {
3457         cused++;
3458       }
3459       else{
3460         cnused++;
3461       }
3462       Double_t x = kr0.GetX();
3463       // Initialization of the polytrack
3464       nfound =0;
3465       nfoundable =0;
3466       polytrack.Reset();
3467       //
3468       Double_t y0= cl->GetY();
3469       Double_t z0= cl->GetZ();
3470       Float_t erry = 0;
3471       Float_t errz = 0;
3472       
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
3475       
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);
3479
3480       sumused=0;
3481       if (cl->IsUsed(10)) sumused++;
3482
3483
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;
3486       //
3487       x = krm.GetX();
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);
3494       }
3495       //
3496       x = krp.GetX();
3497       AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3498       if (cl2) {
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);
3503       }
3504       //
3505       if (sumused>0) continue;
3506       nin0++;
3507       polytrack.UpdateParameters();
3508       // follow polytrack
3509       roadz = 1.2;
3510       roady = 1.2;
3511       //
3512       Double_t yn,zn;
3513       nfoundable = polytrack.GetN();
3514       nfound     = nfoundable; 
3515       //
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;
3525           nfoundable++;
3526           AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3527           if (cln) {
3528             Float_t dist =  TMath::Sqrt(  (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3529             if (dist<maxdist){
3530               /*
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");
3535                 sumused++;
3536                 erry*=2;
3537                 errz*=2;
3538               }
3539               */
3540               erry=0.1;
3541               errz=0.1;
3542               polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3543               nfound++;
3544             }
3545           }
3546         }
3547         if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable))  break;     
3548         polytrack.UpdateParameters();
3549       }           
3550     }
3551     if ( (sumused>3) || (sumused>0.5*nfound))  {
3552       //printf("sumused   %d\n",sumused);
3553       continue;
3554     }
3555     nin1++;
3556     Double_t dy,dz;
3557     polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3558     AliTPCpolyTrack track2;
3559     
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;
3562     nin2++;
3563
3564     if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3565       //
3566       // test seed with and without constrain
3567       for (Int_t constrain=0; constrain<=0;constrain++){
3568         // add polytrack candidate
3569
3570         Double_t x[5], c[15];
3571         Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3572         track2.GetBoundaries(x3,x1);    
3573         x2 = (x1+x3)/2.;
3574         track2.GetFitPoint(x1,y1,z1);
3575         track2.GetFitPoint(x2,y2,z2);
3576         track2.GetFitPoint(x3,y3,z3);
3577         //
3578         //is track pointing to the vertex ?
3579         Double_t x0,y0,z0;
3580         x0=0;
3581         polytrack.GetFitPoint(x0,y0,z0);
3582
3583         if (constrain) {
3584           x2 = x3;
3585           y2 = y3;
3586           z2 = z3;
3587           
3588           x3 = 0;
3589           y3 = 0;
3590           z3 = 0;
3591         }
3592         x[0]=y1;
3593         x[1]=z1;
3594         x[4]=F1(x1,y1,x2,y2,x3,y3);
3595                 
3596         //      if (TMath::Abs(x[4]) >= cuts[0]) continue;  //
3597         x[2]=F2(x1,y1,x2,y2,x3,y3);
3598         
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;
3603
3604         
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;
3608         Double_t sy3=0.02;
3609
3610         if (constrain){
3611           sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3612         }
3613         
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;
3620
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;
3625
3626         
3627         c[0]=sy1;
3628         c[1]=0.;       c[2]=sz1;
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;
3635         
3636         //Int_t row1 = fSectors->GetRowNumber(x1);
3637         Int_t row1 = GetRowNumber(x1);
3638
3639         UInt_t index=0;
3640         //kr0.GetIndex(is);
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);
3646         else
3647           track->SetBConstrain(0);
3648         track->SetLastPoint(row1+fInnerSec->GetNRows());  // first cluster in track position
3649         track->SetFirstPoint(track->GetLastPoint());
3650
3651         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 || 
3652             track->GetNumberOfClusters() < track->GetNFoundable()*0.6 || 
3653             track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3654           //delete track;
3655           seed->Reset();
3656           seed->~AliTPCseed();
3657         }
3658         else {
3659           arr->AddLast(track);
3660           seed = new AliTPCseed;
3661         }
3662         nin3++;
3663       }
3664     }  // if accepted seed
3665   }
3666   if (fDebug>3){
3667     Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3668   }
3669   delete seed;
3670 }
3671
3672
3673 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3674 {
3675   //
3676   //
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
3681   Int_t pp2=0;
3682   Double_t  x0[3],x1[3],x2[3];
3683   for (Int_t i=0;i<3;i++){
3684     x0[i]=-1;
3685     x1[i]=-1;
3686     x2[i]=-1;
3687   }
3688
3689   // find track position at given ratio of the length
3690   Int_t  sec0=0, sec1=0, sec2=0;
3691   Int_t index=-1;
3692   Int_t clindex;
3693   for (Int_t i=0;i<160;i++){
3694     if (track->GetClusterPointer(i)){
3695       index++;
3696       AliTPCTrackerPoint   *trpoint =track->GetTrackPoint(i);
3697       if ( (index<p0) || x0[0]<0 ){
3698         if (trpoint->GetX()>1){
3699           clindex = track->GetClusterIndex2(i);
3700           if (clindex>0){       
3701             x0[0] = trpoint->GetX();
3702             x0[1] = trpoint->GetY();
3703             x0[2] = trpoint->GetZ();
3704             sec0  = ((clindex&0xff000000)>>24)%18;
3705           }
3706         }
3707       }
3708
3709       if ( (index<p1) &&(trpoint->GetX()>1)){
3710         clindex = track->GetClusterIndex2(i);
3711         if (clindex>0){
3712           x1[0] = trpoint->GetX();
3713           x1[1] = trpoint->GetY();
3714           x1[2] = trpoint->GetZ();
3715           sec1  = ((clindex&0xff000000)>>24)%18;
3716         }
3717       }
3718       if ( (index<p2) &&(trpoint->GetX()>1)){
3719         clindex = track->GetClusterIndex2(i);
3720         if (clindex>0){
3721           x2[0] = trpoint->GetX();
3722           x2[1] = trpoint->GetY();
3723           x2[2] = trpoint->GetZ(); 
3724           sec2  = ((clindex&0xff000000)>>24)%18;
3725           pp2 = i;
3726         }
3727       }
3728     }
3729   }
3730   
3731   Double_t alpha, cs,sn, xx2,yy2;
3732   //
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;
3738   x1[0] = xx2;
3739   x1[1] = yy2;
3740   //
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;
3746   x0[0] = xx2;
3747   x0[1] = yy2;
3748   //
3749   //
3750   //
3751   Double_t x[5],c[15];
3752   //
3753   x[0]=x2[1];
3754   x[1]=x2[2];
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;
3761   //  
3762   Double_t sy =0.1,  sz =0.1;
3763   //
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();
3767   //
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;
3774   //
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;
3779   
3780   
3781   c[0]=sy1;
3782   c[1]=0.;       c[2]=sz1;
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;
3789   
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);
3796   //  seed =0;
3797   seed->SetLastPoint(pp2);
3798   seed->SetFirstPoint(pp2);
3799   
3800
3801   return seed;
3802 }
3803
3804
3805 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3806 {
3807   //
3808   //
3809   //reseed using founded clusters 
3810   //
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++;
3815   }
3816   //
3817   Int_t ipos[3];
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
3821   //
3822   //
3823   Double_t  xyz[3][3];
3824   Int_t     row[3],sec[3]={0,0,0};
3825   //
3826   // find track row position at given ratio of the length
3827   Int_t index=-1;
3828   for (Int_t irow=0;irow<160;irow++){    
3829     if (track->GetClusterIndex2(irow)<0) continue;
3830     index++;
3831     for (Int_t ipoint=0;ipoint<3;ipoint++){
3832       if (index<=ipos[ipoint]) row[ipoint] = irow;
3833     }        
3834   }
3835   //
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);
3840     if (cl==0) {
3841       //Error("Bug\n");
3842       //      AliTPCclusterMI * cl = GetClusterMI(clindex);
3843       return 0;
3844     }
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();
3849   }
3850   //
3851   //
3852   // Calculate seed state vector and covariance matrix
3853
3854   Double_t alpha, cs,sn, xx2,yy2;
3855   //
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;
3861   xyz[1][0] = xx2;
3862   xyz[1][1] = yy2;
3863   //
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;
3869   xyz[0][0] = xx2;
3870   xyz[0][1] = yy2;
3871   //
3872   //
3873   //
3874   Double_t x[5],c[15];
3875   //
3876   x[0]=xyz[2][1];
3877   x[1]=xyz[2][2];
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]);
3881   //  
3882   Double_t sy =0.1,  sz =0.1;
3883   //
3884   Double_t sy1=0.2, sz1=0.2;
3885   Double_t sy2=0.2, sz2=0.2;
3886   Double_t sy3=0.2;
3887   //
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;
3894   //
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;
3899   
3900   
3901   c[0]=sy1;
3902   c[1]=0.;       c[2]=sz1;
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;
3909   
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]);  
3914   return seed;
3915 }
3916
3917
3918 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3919 {
3920   //
3921   //
3922   //reseed using founded clusters 
3923   //
3924   Double_t  xyz[3][3];
3925   Int_t     row[3]={0,0,0};
3926   Int_t     sec[3]={0,0,0};
3927   //
3928   // forward direction
3929   if (forward){
3930     for (Int_t irow=r0;irow<160;irow++){
3931       if (track->GetClusterIndex(irow)>0){
3932         row[0] = irow;
3933         break;
3934       }
3935     }
3936     for (Int_t irow=160;irow>r0;irow--){
3937       if (track->GetClusterIndex(irow)>0){
3938         row[2] = irow;
3939         break;
3940       }
3941     }
3942     for (Int_t irow=row[2]-15;irow>row[0];irow--){
3943       if (track->GetClusterIndex(irow)>0){
3944         row[1] = irow;
3945         break;
3946       }
3947     }
3948     //
3949   }
3950   if (!forward){
3951     for (Int_t irow=0;irow<r0;irow++){
3952       if (track->GetClusterIndex(irow)>0){
3953         row[0] = irow;
3954         break;
3955       }
3956     }
3957     for (Int_t irow=r0;irow>0;irow--){
3958       if (track->GetClusterIndex(irow)>0){
3959         row[2] = irow;
3960         break;
3961       }
3962     }    
3963     for (Int_t irow=row[2]-15;irow>row[0];irow--){
3964       if (track->GetClusterIndex(irow)>0){
3965         row[1] = irow;
3966         break;
3967       }
3968     } 
3969   }
3970   //
3971   if ((row[2]-row[0])<20) return 0;
3972   if (row[1]==0) return 0;
3973   //
3974   //
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);
3979     if (cl==0) {
3980       //Error("Bug\n");
3981       //      AliTPCclusterMI * cl = GetClusterMI(clindex);
3982       return 0;
3983     }
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){
3988       //
3989        xyz[ipoint][1]  = point->GetY();
3990        xyz[ipoint][2]  = point->GetZ();
3991     }
3992     else{
3993       xyz[ipoint][1]  = cl->GetY();
3994       xyz[ipoint][2]  = cl->GetZ();
3995     }
3996   }
3997   //
3998   //
3999   //
4000   //
4001   // Calculate seed state vector and covariance matrix
4002
4003   Double_t alpha, cs,sn, xx2,yy2;
4004   //
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;
4010   xyz[1][0] = xx2;
4011   xyz[1][1] = yy2;
4012   //
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;
4018   xyz[0][0] = xx2;
4019   xyz[0][1] = yy2;
4020   //
4021   //
4022   //
4023   Double_t x[5],c[15];
4024   //
4025   x[0]=xyz[2][1];
4026   x[1]=xyz[2][2];
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]);
4030   //  
4031   Double_t sy =0.1,  sz =0.1;
4032   //
4033   Double_t sy1=0.2, sz1=0.2;
4034   Double_t sy2=0.2, sz2=0.2;
4035   Double_t sy3=0.2;
4036   //
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;
4043   //
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;
4048   
4049   
4050   c[0]=sy1;
4051   c[1]=0.;       c[2]=sz1;
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;
4058   
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));
4065   }
4066
4067   return seed;
4068 }
4069
4070
4071
4072 void  AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4073 {
4074   //
4075   //  find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4076   //                      USES MC LABELS
4077   //  Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4078   //
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
4084   //
4085   //
4086   //  Algorithm done in 2 phases - because of CPU consumption
4087   //  it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated                           
4088   //
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)
4094   //                 d.) delta phi
4095   //    2 phase - calculates DCA between tracks  - time consument
4096
4097   //
4098   //    fast cuts 
4099   //
4100   //    General cuts    - for splitted tracks and for curling tracks
4101   //
4102   const Float_t kMaxdPhi      = 0.2;  // maximal distance in phi
4103   //
4104   //    Curling tracks cuts
4105   //
4106   //
4107   //
4108   //
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];
4114   //
4115   //
4116   TStopwatch timer;
4117   timer.Start();
4118   //
4119   // Find track COG in x direction - point with best defined parameters
4120   //
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);
4126     Int_t ncl=0;
4127     xm[i]=0;
4128     Float_t dz[2];
4129     track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4130     dz0[i]=dz[0];
4131     dz1[i]=dz[1];
4132     for (Int_t icl=0; icl<160; icl++){
4133       AliTPCclusterMI * cl =  track->GetClusterPointer(icl);
4134       if (cl) {
4135         xm[i]+=cl->GetX();
4136         ncl++;
4137       }
4138     }
4139     if (ncl>0) xm[i]/=Float_t(ncl);
4140   }  
4141   TTreeSRedirector &cstream = *fDebugStreamer;
4142   //
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);
4151     
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;
4158       //
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);
4164       //
4165       Float_t dfi = fi0-fi1;
4166       //
4167       //
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){
4171         //
4172         // if short tracks with undefined sign 
4173         fi1 =  -TMath::ATan2(yc1,-xc1);
4174         dfi = fi0-fi1;
4175       }
4176       Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4177       
4178       //
4179       // debug stream to tune "fast cuts" 
4180       //
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;
4190       
4191       Float_t sum =0;
4192       Float_t sums=0;
4193       for (Int_t icl=0; icl<160; icl++){
4194         AliTPCclusterMI * cl0 =  track0->GetClusterPointer(icl);
4195         AliTPCclusterMI * cl1 =  track1->GetClusterPointer(icl);
4196         if (cl0&&cl1) {
4197           sum++;
4198           if (cl0==cl1) sums++;
4199         }
4200       }
4201       //
4202       cstream<<"Multi"<<
4203         "iter="<<iter<<
4204         "lab0="<<lab0<<
4205         "lab1="<<lab1<<   
4206         "Tr0.="<<track0<<       // seed0
4207         "Tr1.="<<track1<<       // seed1
4208         "h0.="<<&helixes[i0]<<
4209         "h1.="<<&helixes[i1]<<
4210         //
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