Tracking SA. Bug fix and spped impovements (E. Crescio)
[u/mrichter/AliRoot.git] / ITS / AliITStrackerSA.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2003, 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 /* $Id$ */
17
18 ////////////////////////////////////////////////////
19 //  Stand alone tracker class                     //
20 //  Origin:  Elisabetta Crescio                   //
21 //  e-mail:  crescio@to.infn.it                   //
22 //  tracks are saved as AliITStrackV2 objects     //
23 ////////////////////////////////////////////////////
24
25 #include <stdlib.h>
26
27 #include <TArrayI.h>
28 #include <TBranch.h>
29 #include <TObjArray.h>
30 #include <TTree.h>
31
32 #include "AliESD.h"
33 #include "AliESDVertex.h"
34 #include "AliESDtrack.h"
35 #include "AliITSVertexer.h"
36 #include "AliITSclusterTable.h"
37 #include "AliITSclusterV2.h"
38 #include "AliITSgeom.h"
39 #include "AliITStrackSA.h"
40 #include "AliITStrackerSA.h"
41 #include "AliRun.h"
42
43 ClassImp(AliITStrackerSA)
44
45 //____________________________________________________________________________
46 AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(){
47   // Default constructor
48   Init();
49  
50 }
51 //____________________________________________________________________________
52 AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom):AliITStrackerMI(geom) 
53 {
54   // Standard constructor (Vertex is known and passed to this obj.)
55   Init();
56   fVert = 0;
57   fGeom = geom;
58  
59 }
60
61 //____________________________________________________________________________
62 AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliESDVertex *vert):AliITStrackerMI(geom) 
63 {
64   // Standard constructor (Vertex is known and passed to this obj.)
65   Init();
66   fVert = vert;
67   fGeom = geom;
68  
69 }
70
71 //____________________________________________________________________________
72 AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliITSVertexer *vertexer):AliITStrackerMI(geom) 
73 {
74   // Standard constructor (Vertex is unknown - vertexer is passed to this obj)
75   Init();
76   fVertexer = vertexer;
77   fGeom = geom;
78  
79 }
80
81 //____________________________________________________________________________
82 AliITStrackerSA::AliITStrackerSA(const AliITStrackerSA& tracker):AliITStrackerMI(){
83   // Copy constructor
84   fPhiEstimate = tracker.fPhiEstimate;
85   for(Int_t i=0;i<2;i++){
86     fPoint1[i]=tracker.fPoint1[i];
87     fPoint2[i]=tracker.fPoint2[i];
88     fPoint3[i]=tracker.fPoint3[i];
89     fPointc[i]=tracker.fPointc[i];
90   }
91   fLambdac = tracker.fLambdac;
92   fPhic = tracker.fPhic;
93   fCoef1 = tracker.fCoef1;
94   fCoef2 = tracker.fCoef2;
95   fCoef3 = tracker.fCoef3;
96   fNloop = tracker.fNloop;
97   fPhiWin = tracker.fPhiWin;
98   fLambdaWin = tracker.fLambdaWin;
99   if(tracker.fVertexer && tracker.fVert){
100     fVert = new AliESDVertex(*tracker.fVert);
101   }
102   else {
103     fVert = tracker.fVert;
104   }
105   fVertexer = tracker.fVertexer;
106   fGeom = tracker.fGeom;
107   fListOfTracks = tracker.fListOfTracks;
108   for(Int_t i=0;i<fGeom->GetNlayers();i++){
109     fCluLayer[i] = tracker.fCluLayer[i];
110     fCluCoord[i] = tracker.fCluCoord[i];
111   } 
112 }
113 //______________________________________________________________________
114 AliITStrackerSA& AliITStrackerSA::operator=(const AliITStrackerSA& /*source*/){
115     // Assignment operator. This is a function which is not allowed to be
116     // done.
117     Error("operator=","Assignment operator not allowed\n");
118     return *this; 
119 }
120
121 //____________________________________________________________________________
122 AliITStrackerSA::~AliITStrackerSA(){
123   // destructor
124   // if fVertexer is not null, the AliESDVertex obj. is owned by this class
125   // and is deleted here
126   if(fVertexer){
127     if(fVert)delete fVert;
128   }
129   fVert = 0;
130   fVertexer = 0;
131  
132   if(fPhiWin)delete []fPhiWin;
133   if(fLambdaWin)delete []fLambdaWin;
134   fListOfTracks->Delete();
135   if(fCluLayer){
136     for(Int_t i=0;i<fGeom->GetNlayers();i++){
137       if(fCluLayer[i]){
138         fCluLayer[i]->Delete();
139         delete fCluLayer[i];
140       }
141     }
142     delete [] fCluLayer;
143   }
144   if(fCluCoord){
145     for(Int_t i=0;i<fGeom->GetNlayers();i++){
146       if(fCluCoord[i]){
147         fCluCoord[i]->Delete();
148         delete fCluCoord[i];
149       }
150     }
151     delete [] fCluCoord;
152   }
153   
154   fGeom = 0; 
155 }
156
157 //____________________________________________________________________________
158 void AliITStrackerSA::Init(){
159   //  Reset all data members
160     fPhiEstimate=0;
161     for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
162     fLambdac=0;
163     fPhic=0;
164     fCoef1=0;
165     fCoef2=0;
166     fCoef3=0;
167     fPointc[0]=0;
168     fPointc[1]=0;
169     fVert = 0;
170     fVertexer = 0;
171     fGeom = 0;
172     SetWindowSizes();
173     fITSclusters = 0;
174     SetSixPoints();
175     fListOfTracks=new TObjArray(0,0);
176     fCluLayer = 0;
177     fCluCoord = 0;
178  }
179 //_______________________________________________________________________
180 void AliITStrackerSA::ResetForFinding(){
181   //  Reset data members used in all loops during track finding
182     fPhiEstimate=0;
183     for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
184     fLambdac=0;
185     fPhic=0;
186     fCoef1=0;
187     fCoef2=0;
188     fCoef3=0;
189     fPointc[0]=0;
190     fPointc[1]=0;
191     fListOfTracks->Delete();
192 }
193
194  
195
196 //______________________________________________________________________
197 Int_t AliITStrackerSA::FindTracks(AliESD* event){
198
199 // Track finder using the ESD object
200
201
202   //controllare numero cluster sui layer1 e 2 (morti?)
203   //non trova tracce...controllare..
204
205   if(!fITSclusters){
206     Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
207     return -1;
208   }
209   
210    
211   //Reads event and mark clusters of traks already found, with flag kITSin
212    Int_t nentr=event->GetNumberOfTracks();
213    while (nentr--) {
214      AliESDtrack *track=event->GetTrack(nentr);
215      if (track->GetStatus()&AliESDtrack::kITSin==AliESDtrack::kITSin){
216        UInt_t idx[6];
217        Int_t ncl = track->GetITSclusters(idx);
218        for(Int_t k=0;k<ncl;k++){
219          AliITSclusterV2* cll = (AliITSclusterV2*)GetCluster(idx[k]);
220          cll->SetBit(kSAflag);
221        }
222      }
223    }
224
225    Double_t primaryVertex[3];
226    event->GetVertex()->GetXYZ(primaryVertex);
227    //Creates TClonesArray with clusters for each layer. The clusters already used
228    //by AliITStrackerMI are not considered
229    
230    Int_t nclusters[6]={0,0,0,0,0,0};
231    Int_t dmar[6]={0,0,0,0,0,0};
232    fCluLayer = new TClonesArray*[fGeom->GetNlayers()];
233    fCluCoord = new TClonesArray*[fGeom->GetNlayers()];
234
235    for(Int_t i=0;i<fGeom->GetNlayers();i++){
236      AliITSlayer &layer=fgLayers[i];
237      for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
238        AliITSclusterV2* cls = (AliITSclusterV2*)layer.GetCluster(cli);
239        if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
240        if(cls->GetQ()==0) continue; //fake clusters dead zones
241        nclusters[i]++;
242      }
243      dmar[i]=0;
244      fCluLayer[i] = new TClonesArray("AliITSclusterV2",nclusters[i]);
245      fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
246    }
247
248    
249    Int_t * firstmod = new Int_t[fGeom->GetNlayers()];
250    for(Int_t i=0;i<fGeom->GetNlayers();i++){
251      firstmod[i]=fGeom->GetModuleIndex(i+1,1,1);
252    }
253    
254    for(Int_t ilay=0;ilay<fGeom->GetNlayers();ilay++){
255      TClonesArray &clulay = *fCluLayer[ilay];
256      TClonesArray &clucoo = *fCluCoord[ilay];
257      AliITSlayer &layer=fgLayers[ilay];
258      for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
259        AliITSclusterV2* cls = (AliITSclusterV2*)layer.GetCluster(cli);
260        if(cls->TestBit(kSAflag)==kTRUE) continue;
261        if(cls->GetQ()==0) continue;
262        Double_t phi=0;Double_t lambda=0;
263        Float_t x=0;Float_t y=0;Float_t z=0;
264        Float_t sx=0;Float_t sy=0;Float_t sz=0;
265        Int_t module = cls->GetDetectorIndex()+firstmod[ilay];
266        GetCoorAngles(cls,module,phi,lambda,x,y,z,primaryVertex);
267        GetCoorErrors(cls,module,sx,sy,sz);
268        new (clulay[dmar[ilay]]) AliITSclusterV2(*cls);
269        new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
270        dmar[ilay]++;
271      }
272    }
273    
274  
275   //Get primary vertex
276
277    Int_t ntrack=0;
278    //loop on the different windows
279    for(Int_t nloop=0;nloop<fNloop;nloop++){
280      for(Int_t ncl=0;ncl<fCluLayer[0]->GetEntries();ncl++){ //loop starting from layer 0
281        
282        ResetForFinding();
283        Int_t pflag=0;
284        
285        AliITSclusterV2* cl = (AliITSclusterV2*)fCluLayer[0]->At(ncl);
286
287        if(!cl) continue;
288
289
290        if (cl->GetQ()<=0) continue;
291        
292        AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(0,ncl); 
293        fPhic = arr->GetPhi();
294        fLambdac = arr->GetLambda();
295        if (TMath::Abs(fLambdac)>0.26*TMath::Pi()) continue;
296        fPhiEstimate = fPhic;
297        AliITStrackSA* trs = new AliITStrackSA(); 
298        fPoint1[0]=primaryVertex[0];
299        fPoint1[1]=primaryVertex[1];
300
301
302        fPoint2[0]=arr->GetX();
303        fPoint2[1]=arr->GetY();
304        Int_t * nn = new Int_t[fGeom->GetNlayers()];//counter for clusters on each layer
305        for(Int_t i=0;i<fGeom->GetNlayers();i++){ nn[i]=0;}
306        nn[0] = SearchClusters(0,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
307        
308        nn[1] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
309       if(nn[1]>0){
310         pflag=1;
311         fPoint3[0] = fPointc[0];
312         fPoint3[1] = fPointc[1];
313       }
314       nn[2] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
315       if(nn[1]==0 && nn[2]==0) pflag=0;
316       if(nn[2]!=0 && nn[1]!=0){ pflag=1; UpdatePoints();}
317       if(nn[2]!=0 && nn[1]==0){
318         pflag=1;
319         fPoint3[0]=fPointc[0];
320         fPoint3[1]=fPointc[1];
321       }
322
323       nn[3] = SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
324       pflag=1;
325       if(nn[3]!=0) UpdatePoints();
326       nn[4] = SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag); 
327       pflag=1;
328       if(nn[4]!=0) UpdatePoints();
329       nn[5] = SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag); 
330           
331
332       Int_t layOK=0;
333       Int_t numberofpoints;
334       if(fSixPoints) numberofpoints=6;  //check of the candidate track
335       else numberofpoints=5;           //if track is good (with the required number        
336       for(Int_t nnp=0;nnp<fGeom->GetNlayers();nnp++){    //of points) it is written on file
337         if(nn[nnp]!=0) layOK+=1;
338       }
339       if(layOK>=numberofpoints){
340
341         AliITStrackV2* tr2 = FitTrack(trs,primaryVertex);
342         
343         if(tr2==0) continue;
344
345         AliESDtrack outtrack;
346         outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
347         event->AddTrack(&outtrack);
348         ntrack++;
349         
350       }
351         
352
353       delete trs;
354       delete[] nn;
355       
356      }//end loop on clusters of layer1
357     
358      //end loop2
359    }
360    
361    //if 5/6 points are required, second loop starting 
362    //from second layer, to find tracks with point of 
363    //layer 1 missing
364    
365    if(!fSixPoints){
366      //   counter for clusters on each layer  
367      Int_t * nn = new Int_t[fGeom->GetNlayers()-1];      
368      for(Int_t nloop=0;nloop<fNloop;nloop++){
369        Int_t ncl2=fCluLayer[1]->GetEntries();
370        while(ncl2--){ //loop starting from layer 2
371          ResetForFinding();
372          Int_t pflag=0;
373          AliITSclusterV2* cl = (AliITSclusterV2*)fCluLayer[1]->At(ncl2);
374          
375          if(!cl) continue;
376          AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(1,ncl2);
377          fPhic = arr->GetPhi();
378          fLambdac = arr->GetLambda();
379          fPhiEstimate = fPhic;
380          
381          AliITStrackSA* trs = new AliITStrackSA(); 
382          fPoint1[0]=primaryVertex[0];
383          fPoint1[1]=primaryVertex[1];
384          
385          fPoint2[0]=arr->GetX();
386          fPoint2[1]=arr->GetY();
387          for(Int_t kk=0;kk<fGeom->GetNlayers()-1;kk++)nn[kk] = 0;
388          nn[0] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],
389                                    trs,primaryVertex[2],pflag);
390          nn[1] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],
391                                    trs,primaryVertex[2],pflag);
392          if(nn[1]!=0){
393            pflag=1;
394            fPoint3[0]=fPointc[0];
395            fPoint3[1]=fPointc[1];
396          }
397          nn[2]= SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],
398                                    trs,primaryVertex[2],pflag);
399          if(nn[2]!=0){
400            pflag=1;
401            UpdatePoints();
402          }
403          nn[3]= SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],
404                                    trs,primaryVertex[2],pflag);
405          if(nn[3]!=0){
406            pflag=1;
407            UpdatePoints();
408          }
409          nn[4]=SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],
410                                    trs,primaryVertex[2],pflag);
411
412          Int_t fl=0;
413          for(Int_t nnp=0;nnp<fGeom->GetNlayers()-1;nnp++){
414            if(nn[nnp]!=0) fl+=1;
415          }
416          if(fl>=5){  // 5/6       
417            AliITStrackV2* tr2 = FitTrack(trs,primaryVertex);
418            if(tr2==0){
419              continue;
420            }
421            
422            AliESDtrack outtrack;
423            outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
424            event->AddTrack(&outtrack);
425            ntrack++;
426
427          }   
428         
429        
430         delete trs;
431       }//end loop on clusters of layer2
432      }
433    
434     delete [] nn;
435   }  //end opt="5/6"  
436    
437   delete [] firstmod;
438   Info("FindTracks","Number of found tracks: %d",event->GetNumberOfTracks());
439   return 0;
440
441 }
442  
443
444
445
446
447 //________________________________________________________________________
448
449 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex){
450   //fit of the found track
451
452   
453   Int_t * firstmod = new Int_t[fGeom->GetNlayers()];
454   for(Int_t i=0;i<fGeom->GetNlayers();i++){
455     firstmod[i]=fGeom->GetModuleIndex(i+1,1,1);
456   }  
457
458   Int_t nclusters = tr->GetNumberOfClustersSA();
459   TObjArray** listlayer = new TObjArray*[fGeom->GetNlayers()];
460   for(Int_t i=0;i<fGeom->GetNlayers();i++){
461     listlayer[i] = new TObjArray(0,0);
462   }
463
464   TArrayI clind0(20);
465   TArrayI clind1(20);
466   TArrayI clind2(20);
467   TArrayI clind3(20);
468   TArrayI clind4(20);
469   TArrayI clind5(20);
470
471   TArrayI mark0(20);
472   TArrayI mark1(20);
473   TArrayI mark2(20);
474   TArrayI mark3(20);
475   TArrayI mark4(20);
476   TArrayI mark5(20);
477
478
479   Int_t * nnn = new Int_t[fGeom->GetNlayers()];
480   Int_t * kkk = new Int_t[fGeom->GetNlayers()];
481   for(Int_t i=0;i<fGeom->GetNlayers();i++) {nnn[i]=0;kkk[i]=0;}
482   
483   for(Int_t ncl=0;ncl<nclusters;ncl++){
484     Int_t index = tr->GetClusterIndexSA(ncl); 
485     AliITSclusterV2* cl = (AliITSclusterV2*)GetCluster(index);
486     if(cl->TestBit(kSAflag)==kTRUE) cl->ResetBit(kSAflag);
487     Int_t lay = (index & 0xf0000000) >> 28;
488     if(lay==0) { listlayer[0]->AddLast(cl); clind0[nnn[0]]=index;nnn[0]++;}
489     if(lay==1) { listlayer[1]->AddLast(cl); clind1[nnn[1]]=index;nnn[1]++;}
490     if(lay==2) { listlayer[2]->AddLast(cl); clind2[nnn[2]]=index;nnn[2]++;}
491     if(lay==3) { listlayer[3]->AddLast(cl); clind3[nnn[3]]=index;nnn[3]++;}
492     if(lay==4) { listlayer[4]->AddLast(cl); clind4[nnn[4]]=index;nnn[4]++;}
493     if(lay==5) { listlayer[5]->AddLast(cl); clind5[nnn[5]]=index;nnn[5]++;}    
494   }
495   delete [] nnn;
496
497   for(Int_t nlay=0;nlay<fGeom->GetNlayers();nlay++){
498     for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
499       Int_t mark = tr->GetClusterMark(nlay,ncl);
500       if(nlay==0) { mark0[kkk[0]]=mark;kkk[0]++;}
501       if(nlay==1) { mark1[kkk[1]]=mark;kkk[1]++;}
502       if(nlay==2) { mark2[kkk[2]]=mark;kkk[2]++;}
503       if(nlay==3) { mark3[kkk[3]]=mark;kkk[3]++;}
504       if(nlay==4) { mark4[kkk[4]]=mark;kkk[4]++;}
505       if(nlay==5) { mark5[kkk[5]]=mark;kkk[5]++;}
506
507     }
508   }
509
510   delete [] kkk;
511
512  
513   Int_t * end = new Int_t[fGeom->GetNlayers()];
514   for(Int_t i=0;i<fGeom->GetNlayers();i++){
515     if(listlayer[i]->GetEntries()==0) end[i]=1;
516     else end[i]=listlayer[i]->GetEntries();
517   }
518
519   TClonesArray* listSA = new TClonesArray("AliITStrackSA");
520   TClonesArray &tri = *listSA;
521   Int_t nlist=0;
522
523   if(end[0]==0) end[0]=1; //for tracks with cluster on layer 0 missing
524   for(Int_t l1=0;l1<end[0];l1++){//loop on layer 1
525     AliITSclusterV2* cl0 = (AliITSclusterV2*)listlayer[0]->At(l1); 
526     Double_t x1,y1,z1,sx1,sy1,sz1;
527     Double_t x2,y2,z2,sx2,sy2,sz2;
528     AliITSclusterV2* p1=0;
529     AliITSclusterV2* p2=0;
530     Int_t index1=clind0[l1];
531     Int_t index2=0;
532     Int_t mrk1 = mark0[l1];
533     Int_t mrk2 = 0;
534     Int_t lay1=0;
535     Int_t lay2=1;
536     Int_t module1=-1;
537     for(Int_t l2=0;l2<end[1];l2++){//loop on layer 2
538       AliITSclusterV2* cl1 = (AliITSclusterV2*)listlayer[1]->At(l2); 
539       index2=clind1[l2];
540       mrk2 = mark1[l2];
541       for(Int_t l3=0;l3<end[2];l3++){  //loop on layer 3
542         AliITSclusterV2* cl2 = (AliITSclusterV2*)listlayer[2]->At(l3);
543
544         if(cl0==0 && cl1!=0) {
545           p2 = cl2;index1=clind2[l3];mrk1=mark2[l3];lay1=2;
546           p1=cl1;
547           module1 = p1->GetDetectorIndex()+firstmod[1]; 
548         }
549         if(cl0!=0 && cl1==0){
550           p1=cl0;
551           p2=cl2;index2=clind2[l3];mrk2=mark2[l3];lay2=2;
552           module1 = p1->GetDetectorIndex()+firstmod[0];
553         }
554         if(cl0!=0 && cl1!=0){
555           p1=cl0;
556           p2=cl1;
557           module1 = p1->GetDetectorIndex()+firstmod[0];
558         }
559         
560         Int_t cln1=mrk1;
561         Int_t cln2=mrk2;
562         AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay1,cln1);
563         AliITSclusterTable* arr1= (AliITSclusterTable*)GetClusterCoord(lay2,cln2);
564         x1 = arr->GetX();
565         x2 = arr1->GetX();
566         y1 = arr->GetY();
567         y2 = arr1->GetY();
568         z1 = arr->GetZ();
569         z2 = arr1->GetZ();
570         sx1 = arr->GetSx();
571         sx2 = arr1->GetSx();
572         sy1 = arr->GetSy();
573         sy2 = arr1->GetSy();
574         sz1 = arr->GetSz();
575         sz2 = arr1->GetSz();
576         
577         Int_t layer,ladder,detector;
578         fGeom->GetModuleId(module1,layer,ladder,detector);
579         Float_t yclu1 = p1->GetY();
580         Float_t zclu1 = p1->GetZ();
581         Double_t cv=Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
582         
583         Double_t tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
584         Double_t phi2 = TMath::ATan2((y2-y1),(x2-x1));
585
586         for(Int_t l4=0;l4<end[3];l4++){ //loop on layer 4   
587           AliITSclusterV2* cl3 = (AliITSclusterV2*)listlayer[3]->At(l4);
588           for(Int_t l5=0;l5<end[4];l5++){ //loop on layer 5
589             AliITSclusterV2* cl4 = (AliITSclusterV2*)listlayer[4]->At(l5);
590             for(Int_t l6=0;l6<end[5];l6++){ //loop on layer 6  
591               AliITSclusterV2* cl5 = (AliITSclusterV2*)listlayer[5]->At(l6);
592               AliITStrackSA* trac = new AliITStrackSA(fGeom,layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
593                               
594               if(cl5!=0) {
595                 trac->AddClusterV2(5,(clind5[l6] & 0x0fffffff)>>0);
596                 trac->AddClusterMark(5,mark5[l6]);
597               }
598               if(cl4!=0){
599                 trac->AddClusterV2(4,(clind4[l5] & 0x0fffffff)>>0);
600                 trac->AddClusterMark(4,mark4[l5]);
601               }
602               if(cl3!=0){
603                 trac->AddClusterV2(3,(clind3[l4] & 0x0fffffff)>>0);
604                 trac->AddClusterMark(3,mark3[l4]);
605               }
606               if(cl2!=0){
607                 trac->AddClusterV2(2,(clind2[l3] & 0x0fffffff)>>0);
608                 trac->AddClusterMark(2,mark2[l3]);
609               }
610               if(cl1!=0){
611                 trac->AddClusterV2(1,(clind1[l2] & 0x0fffffff)>>0);
612                 trac->AddClusterMark(1,mark1[l2]);
613               }
614               if(cl0!=0){
615                 trac->AddClusterV2(0,(clind0[l1] & 0x0fffffff)>>0);
616                 trac->AddClusterMark(0,mark0[l1]);
617               }
618               //fit with Kalman filter using AliITStrackerMI::RefitAt()
619           
620
621               AliITStrackMI* ot = new AliITStrackSA(*trac);
622               
623               ot->ResetCovariance();
624               ot->ResetClusters();
625               
626               if(RefitAt(49.,ot,trac)){ //fit from layer 1 to layer 6
627                 AliITStrackMI *otrack2 = new AliITStrackMI(*ot);
628                 otrack2->ResetCovariance(); 
629                 otrack2->ResetClusters();
630                 //fit from layer 6 to layer 1
631                 if(RefitAt(3.7,otrack2,ot)) {
632                   fListOfTracks->AddLast(otrack2);
633                   new (tri[nlist]) AliITStrackSA(*trac);
634                   nlist++;
635                 } else {
636                   delete otrack2;
637                 }
638                               
639               }       
640           
641               delete ot;
642               delete trac;
643             }//end loop layer 6
644           }//end loop layer 5
645         }//end loop layer 4
646         
647       }//end loop layer 3
648     }//end loop layer 2 
649   }//end loop layer 1
650
651   delete [] end;
652
653  
654
655   Int_t dim=fListOfTracks->GetEntries();
656   if(dim==0){
657     for(Int_t i=0;i<fGeom->GetNlayers();i++){
658       delete listlayer[i];
659     }
660     delete [] listlayer;
661     delete [] firstmod;
662     return 0;
663   }
664
665   Int_t lowchi2 = FindTrackLowChiSquare(fListOfTracks,dim);
666   AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
667   AliITStrackSA* trsa = (AliITStrackSA*)listSA->At(lowchi2);
668  
669   if(otrack==0) {
670     for(Int_t i=0;i<fGeom->GetNlayers();i++){
671       delete listlayer[i];
672     }
673     delete [] listlayer; 
674     delete [] firstmod;
675     return 0;
676   }
677   Int_t * indexc = new Int_t[fGeom->GetNlayers()];
678   for(Int_t i=0;i<fGeom->GetNlayers();i++) indexc[i]=0;
679   for(Int_t nind=0;nind<otrack->GetNumberOfClusters();nind++){
680     indexc[nind] = otrack->GetClusterIndex(nind);
681   }      
682   AliITSclusterV2* cl0 = (AliITSclusterV2*)GetCluster(indexc[0]);
683   AliITSclusterV2* cl1 = (AliITSclusterV2*)GetCluster(indexc[1]);     
684   AliITSclusterV2* cl2 = (AliITSclusterV2*)GetCluster(indexc[2]);     
685   AliITSclusterV2* cl3 = (AliITSclusterV2*)GetCluster(indexc[3]);
686   AliITSclusterV2* cl4 = (AliITSclusterV2*)GetCluster(indexc[4]);
687   Int_t labl[3]={-1,-1,-1};
688   if(otrack->GetNumberOfClusters()==fGeom->GetNlayers()){
689     AliITSclusterV2* cl5 = (AliITSclusterV2*)GetCluster(indexc[5]);
690     labl[0]=cl5->GetLabel(0);
691     labl[1]=cl5->GetLabel(1);
692     labl[2]=cl5->GetLabel(2);
693   }
694   delete [] indexc;
695   if(otrack->GetNumberOfClusters()==(fGeom->GetNlayers()-1)){
696     labl[0]=-1;
697     labl[1]=-1;
698     labl[2]=-1;
699   }
700   Int_t numberofpoints;
701   if(fSixPoints) numberofpoints=6;
702   else numberofpoints=5;
703   CookLabel(otrack,0.); //MI change - to see fake ratio
704   Int_t label =  Label(cl0->GetLabel(0),cl1->GetLabel(0), 
705                        cl2->GetLabel(0),cl3->GetLabel(0),
706                        cl4->GetLabel(0),labl[0],
707                        cl0->GetLabel(1),cl1->GetLabel(1),
708                        cl2->GetLabel(1),cl3->GetLabel(1),
709                        cl4->GetLabel(1),labl[1],
710                        cl0->GetLabel(2),cl1->GetLabel(2),
711                        cl2->GetLabel(2),cl3->GetLabel(2),
712                        cl4->GetLabel(2),labl[2],numberofpoints);
713   
714   otrack->SetLabel(label);  
715   //remove clusters of found track
716
717   for(Int_t nlay=0;nlay<fGeom->GetNlayers();nlay++){
718     for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
719       Int_t index = trsa->GetClusterMark(nlay,cln);
720       fCluLayer[nlay]->RemoveAt(index);
721       RemoveClusterCoord(nlay,index);
722       fCluLayer[nlay]->Compress();
723     }    
724   }
725   listSA->Delete();
726   delete listSA;
727
728   for(Int_t i=0;i<fGeom->GetNlayers();i++){
729     delete listlayer[i];
730   }
731   delete [] listlayer; 
732   delete [] firstmod;
733   return otrack;
734
735 }
736
737
738
739 //_______________________________________________________
740 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
741   //function used to to find the clusters associated to the track
742   Int_t nc=0;
743   AliITSlayer &lay = fgLayers[layer];
744   Double_t r=lay.GetR();
745    if(pflag==1){      
746     Float_t cx1,cx2,cy1,cy2;
747     FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
748     if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
749        return 0;
750     Double_t fi1=TMath::ATan2(cy1,cx1);
751     Double_t fi2=TMath::ATan2(cy2,cx2);
752     fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
753   }
754
755  
756   Int_t ncl = fCluLayer[layer]->GetEntries();
757   for (Int_t index=0; index<ncl; index++) {
758     AliITSclusterV2 *c = (AliITSclusterV2*)fCluLayer[layer]->At(index);
759     if (!c) continue;
760     if (c->GetQ()<=0) continue;
761     
762      AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
763      Double_t phi = arr->GetPhi();
764      if (TMath::Abs(phi-fPhiEstimate)>phiwindow) continue;
765
766      Double_t lambda = arr->GetLambda();
767      if (TMath::Abs(lambda-fLambdac)>lambdawindow) continue;
768
769      if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
770      if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
771      Int_t orind = arr->GetOrInd();
772      trs->AddClusterSA(layer,orind);
773      trs->AddClusterMark(layer,index);
774        
775      nc++;
776      fLambdac=lambda;
777      fPhiEstimate=phi;
778
779      fPointc[0]=arr->GetX();
780      fPointc[1]=arr->GetY();
781
782   }
783   return nc;
784 }
785
786 //________________________________________________________________
787 void AliITStrackerSA::UpdatePoints(){
788   //update of points for the estimation of the curvature  
789
790   fPoint2[0]=fPoint3[0];
791   fPoint2[1]=fPoint3[1];
792   fPoint3[0]=fPointc[0];
793   fPoint3[1]=fPointc[1];
794
795   
796 }
797
798 //___________________________________________________________________
799 Int_t AliITStrackerSA::FindEquation(Float_t x1, Float_t y1, Float_t x2, Float_t y2, Float_t x3, Float_t y3,Float_t& a, Float_t& b, Float_t& c){
800
801    //given (x,y) of three recpoints (in global coordinates) 
802    //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
803
804    Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
805    if(den==0) return 0;
806    a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
807    b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
808    c = -x1*x1-y1*y1-a*x1-b*y1;
809    return 1;
810  }
811 //__________________________________________________________________________
812  Int_t AliITStrackerSA::FindIntersection(Float_t a1, Float_t b1, Float_t c1, Float_t c2,Float_t& x1,Float_t& y1, Float_t& x2, Float_t& y2){
813  
814  //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
815  //c2 is -rlayer*rlayer
816
817   if(a1==0) return 0;
818  Double_t m = c2-c1; 
819  Double_t aA = (b1*b1)/(a1*a1)+1;
820  Double_t bB = (-2*m*b1/(a1*a1));
821  Double_t cC = c2+(m*m)/(a1*a1);
822  Double_t dD = bB*bB-4*aA*cC;
823  if(dD<0) return 0;
824  
825  y1 = (-bB+TMath::Sqrt(dD))/(2*aA); 
826  y2 = (-bB-TMath::Sqrt(dD))/(2*aA); 
827  x1 = (c2-c1-b1*y1)/a1;
828  x2 = (c2-c1-b1*y2)/a1;
829
830  return 1; 
831 }
832 //____________________________________________________________________
833 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t 
834 x2,Double_t y2,Double_t x3,Double_t y3){
835
836   //calculates the curvature of track  
837   Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
838   if(den==0) return 0;
839   Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
840   Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
841   Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
842   Double_t xc=-a/2.;
843
844   if((a*a+b*b-4*c)<0) return 0;
845   Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
846   if(rad==0) return 0;
847   
848   if((x1>0 && y1>0 && x1<xc)) rad*=-1;
849   if((x1<0 && y1>0 && x1<xc)) rad*=-1;
850   //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
851   // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
852   
853   return 1/rad;
854  
855 }
856
857
858 //____________________________________________________________________
859 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
860
861   //Returns the point closest to pp
862
863   Double_t diff1 = p1-pp;
864   Double_t diff2 = p2-pp;
865   
866   if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
867   else fPhiEstimate=p2;  
868   return fPhiEstimate;
869   
870 }
871
872
873 //_________________________________________________________________
874 Int_t AliITStrackerSA::FindTrackLowChiSquare(TObjArray* tracklist, Int_t dim) const {
875   // returns track with lowes chi square  
876   if(dim==1){
877     //AliITStrackV2* trk = (AliITStrackV2*)tracklist->At(0);
878     //return trk;
879     return 0;
880   }
881   //if(dim==0) return 0;
882   Double_t * chi2 = new Double_t[dim];
883   Int_t * index = new Int_t[dim];
884   for(Int_t i=0;i<dim;i++){
885     AliITStrackV2* trk = (AliITStrackV2*)tracklist->At(i);
886     chi2[i]=trk->GetChi2();
887     index[i]=i;
888   }
889
890   Int_t w=0;Double_t value;
891   Int_t lp;
892   while(w<dim){
893     for(Int_t j=w+1;j<dim;j++){
894       if(chi2[w]<chi2[j]){
895         value=chi2[w];
896         chi2[w]=chi2[j];
897         chi2[j]=value;
898         lp=index[w];
899         index[w]=index[j];
900         index[j]=lp;
901       }
902     }
903     w++;
904   }
905
906   delete [] chi2;
907   delete [] index;
908   return index[dim-1];
909 }
910
911 //__________________________________________________________
912 Int_t AliITStrackerSA::FindLabel(Int_t l1, Int_t l2, Int_t l3, Int_t l4, Int_t l5, Int_t l6){
913
914   //function used to determine the track label
915   
916   Int_t lb[6] = {l1,l2,l3,l4,l5,l6};
917   Int_t aa[6]={1,1,1,1,1,1};
918   Int_t ff=0; 
919   Int_t ll=0;
920   Int_t k=0;Int_t w=0;Int_t num=6;
921   if(lb[5]==-1) num=5;
922   
923   while(k<num){
924   
925     for(Int_t i=k+1;i<num;i++){
926     
927       if(lb[k]==lb[i] && aa[k]!=0){
928       
929         aa[k]+=1;
930         aa[i]=0;
931       }
932     }
933   k++;
934   }
935
936   while(w<num){
937   
938     for(Int_t j=0;j<6;j++){
939       if(aa[w]<aa[j]){
940       ff=aa[w];
941       aa[w]=aa[j];
942       aa[j]=ff;
943       ll=lb[w];
944       lb[w]=lb[j];
945       lb[j]=ll;
946      }
947     }
948   w++;
949   }
950   if(num==6)  return lb[5];
951   else return lb[4];
952 }
953
954 //_____________________________________________________________________________
955 Int_t AliITStrackerSA::Label(Int_t gl1, Int_t gl2, Int_t gl3, Int_t gl4, Int_t gl5, Int_t gl6,Int_t gl7, Int_t gl8, Int_t gl9, Int_t gl10,Int_t gl11,
956 Int_t gl12, Int_t gl13, Int_t gl14,Int_t gl15, Int_t gl16, Int_t gl17, Int_t gl18, Int_t numberofpoints){
957
958  
959   //function used to assign label to the found track. If track is fake, the label is negative
960
961   Int_t lb0[6] = {gl1,gl2,gl3,gl4,gl5,gl6};
962   Int_t lb1[6] = {gl7,gl8,gl9,gl10,gl11,gl12};
963   Int_t lb2[6] = {gl13,gl14,gl15,gl16,gl17,gl18};
964   Int_t ll=FindLabel(lb0[0],lb0[1],lb0[2],lb0[3],lb0[4],lb0[5]);
965   Int_t lflag=0;Int_t num=6;
966   if(lb0[5]==-1 && lb1[5]==-1 && lb2[5]==-1) num=5;
967
968   for(Int_t i=0;i<num;i++){
969     if(lb0[i]==ll || lb1[i]==ll || lb2[i]==ll) lflag+=1;
970   }
971
972   if(lflag>=numberofpoints) return ll;
973   else return -ll;
974
975   
976 }
977
978 //_____________________________________________________________________________
979 void AliITStrackerSA::SetWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
980   // Set sizes of the phi and lambda windows used for track finding
981   fNloop = n;
982   if(phi){ // user defined values
983     fPhiWin = new Double_t[fNloop];
984     fLambdaWin = new Double_t[fNloop];
985     for(Int_t k=0;k<fNloop;k++){
986       fPhiWin[k]=phi[k];
987       fLambdaWin[k]=lam[k];
988     }
989   }
990   else {  // default values
991             
992     Double_t phid[33]   = {0.002,0.003,0.004,0.0045,0.0047,
993                            0.005,0.0053,0.0055,
994                            0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
995                            0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
996                            0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
997     Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
998                             0.005,0.005,0.006,
999                             0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1000                             0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1001                             0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1002     
1003     if(fNloop!=33){
1004       fNloop = 33;
1005     }
1006     
1007     
1008     fPhiWin = new Double_t[fNloop];
1009     fLambdaWin = new Double_t[fNloop];
1010    
1011     for(Int_t k=0;k<fNloop;k++){
1012       fPhiWin[k]=phid[k];
1013       fLambdaWin[k]=lambdad[k];
1014     }
1015   
1016   }
1017
1018 }
1019 //_______________________________________________________________________
1020 void AliITStrackerSA::GetCoorAngles(AliITSclusterV2* cl,Int_t module,Double_t &phi,Double_t &lambda, Float_t &x, Float_t &y,Float_t &z,Double_t* vertex){
1021   //Returns values of phi (azimuthal) and lambda angles for a given cluster
1022   
1023   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1024   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1025   Float_t tx,ty,tz;  fGeom->GetTrans(lay,lad,det,tx,ty,tz);     
1026
1027   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1028   Double_t phi1=TMath::Pi()/2+alpha;
1029   if (lay==1) phi1+=TMath::Pi();
1030
1031   Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1032   Float_t r=tx*cp+ty*sp;
1033
1034   x= r*cp - cl->GetY()*sp;
1035   y= r*sp + cl->GetY()*cp;
1036   z=cl->GetZ();
1037   
1038   phi=TMath::ATan2(y,x);
1039   lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1040 }
1041
1042 //________________________________________________________________________
1043 void AliITStrackerSA::GetCoorErrors(AliITSclusterV2* cl, Int_t module,Float_t &sx,Float_t &sy, Float_t &sz){
1044
1045   //returns x,y,z of cluster in global coordinates
1046
1047   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1048   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1049  
1050   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1051   Double_t phi=TMath::Pi()/2+alpha;
1052   if (lay==1) phi+=TMath::Pi();
1053
1054   Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1055
1056   sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1057   sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1058   sz = TMath::Sqrt(cl->GetSigmaZ2());
1059
1060 }