warnings with gcc 4.3.2
[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 "AliESDEvent.h"
33 #include "AliESDVertex.h"
34 #include "AliESDtrack.h"
35 #include "AliITSVertexer.h"
36 #include "AliITSclusterTable.h"
37 #include "AliITSRecPoint.h"
38 #include "AliITSgeomTGeo.h"
39 #include "AliITStrackSA.h"
40 #include "AliITStrackerSA.h"
41 #include "AliITSReconstructor.h"
42 #include "AliRun.h"
43
44 ClassImp(AliITStrackerSA)
45
46 //____________________________________________________________________________
47 AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(),
48 fPhiEstimate(0),
49 fITSStandAlone(0),
50 fLambdac(0),
51 fPhic(0),
52 fCoef1(0),
53 fCoef2(0),
54 fCoef3(0),
55 fNloop(0),
56 fPhiWin(0),
57 fLambdaWin(0),
58 fVert(0),
59 fVertexer(0),
60 fListOfTracks(0),
61 fListOfSATracks(0),
62 fITSclusters(0),
63 fSixPoints(0),
64 fOuterStartLayer(0),
65 fCluLayer(0),
66 fCluCoord(0){
67   // Default constructor
68   Init();
69  
70 }
71 //____________________________________________________________________________
72 AliITStrackerSA::AliITStrackerSA(const Char_t *geom):AliITStrackerMI(0),
73 fPhiEstimate(0),
74 fITSStandAlone(0),
75 fLambdac(0),
76 fPhic(0),
77 fCoef1(0),
78 fCoef2(0),
79 fCoef3(0),
80 fNloop(0),
81 fPhiWin(0),
82 fLambdaWin(0),
83 fVert(0),
84 fVertexer(0),
85 fListOfTracks(0),
86 fListOfSATracks(0),
87 fITSclusters(0),
88 fSixPoints(0),
89 fOuterStartLayer(0),
90 fCluLayer(0),
91 fCluCoord(0) 
92 {
93   // Standard constructor (Vertex is known and passed to this obj.)
94   if (geom) {
95     AliWarning("\"geom\" is actually a dummy argument !");
96   }
97
98   Init();
99   fVert = 0;
100  
101 }
102
103 //____________________________________________________________________________
104 AliITStrackerSA::AliITStrackerSA(const Char_t *geom, AliESDVertex *vert):AliITStrackerMI(0),
105 fPhiEstimate(0),
106 fITSStandAlone(0),
107 fLambdac(0),
108 fPhic(0),
109 fCoef1(0),
110 fCoef2(0),
111 fCoef3(0),
112 fNloop(0),
113 fPhiWin(0),
114 fLambdaWin(0),
115 fVert(vert),
116 fVertexer(0),
117 fListOfTracks(0),
118 fListOfSATracks(0),
119 fITSclusters(0),
120 fSixPoints(0),
121 fOuterStartLayer(0),
122 fCluLayer(0),
123 fCluCoord(0)
124 {
125   // Standard constructor (Vertex is known and passed to this obj.)
126   if (geom) {
127     AliWarning("\"geom\" is actually a dummy argument !");
128   }
129   Init();
130  
131 }
132
133 //____________________________________________________________________________
134 AliITStrackerSA::AliITStrackerSA(const Char_t *geom, AliITSVertexer *vertexer):AliITStrackerMI(0),
135 fPhiEstimate(0),
136 fITSStandAlone(0),
137 fLambdac(0),
138 fPhic(0),
139 fCoef1(0),
140 fCoef2(0),
141 fCoef3(0),
142 fNloop(0),
143 fPhiWin(0),
144 fLambdaWin(0),
145 fVert(),
146 fVertexer(vertexer),
147 fListOfTracks(0),
148 fListOfSATracks(0),
149 fITSclusters(0),
150 fSixPoints(0),
151 fOuterStartLayer(0),
152 fCluLayer(0),
153 fCluCoord(0)
154 {
155   // Standard constructor (Vertex is unknown - vertexer is passed to this obj)
156   if (geom) {
157     AliWarning("\"geom\" is actually a dummy argument !");
158   }
159   Init();
160   fVertexer = vertexer;
161  
162 }
163
164 //____________________________________________________________________________
165 AliITStrackerSA::AliITStrackerSA(const AliITStrackerSA& tracker):AliITStrackerMI(),
166 fPhiEstimate(tracker.fPhiEstimate),
167 fITSStandAlone(tracker.fITSStandAlone),
168 fLambdac(tracker.fLambdac),
169 fPhic(tracker.fPhic),
170 fCoef1(tracker.fCoef1),
171 fCoef2(tracker.fCoef2),
172 fCoef3(tracker.fCoef3),
173 fNloop(tracker.fNloop),
174 fPhiWin(tracker.fPhiWin),
175 fLambdaWin(tracker.fLambdaWin),
176 fVert(tracker.fVert),
177 fVertexer(tracker.fVertexer),
178 fListOfTracks(tracker.fListOfTracks),
179 fListOfSATracks(tracker.fListOfSATracks),
180 fITSclusters(tracker.fITSclusters),
181 fSixPoints(tracker.fSixPoints),
182 fOuterStartLayer(tracker.fOuterStartLayer),
183 fCluLayer(tracker.fCluLayer),
184 fCluCoord(tracker.fCluCoord) {
185   // Copy constructor
186   for(Int_t i=0;i<2;i++){
187     fPoint1[i]=tracker.fPoint1[i];
188     fPoint2[i]=tracker.fPoint2[i];
189     fPoint3[i]=tracker.fPoint3[i];
190     fPointc[i]=tracker.fPointc[i];
191   }
192   if(tracker.fVertexer && tracker.fVert){
193     fVert = new AliESDVertex(*tracker.fVert);
194   }
195   else {
196     fVert = tracker.fVert;
197   }
198   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
199     fCluLayer[i] = tracker.fCluLayer[i];
200     fCluCoord[i] = tracker.fCluCoord[i];
201   } 
202 }
203 //______________________________________________________________________
204 AliITStrackerSA& AliITStrackerSA::operator=(const AliITStrackerSA& source){
205     // Assignment operator. 
206   this->~AliITStrackerSA();
207   new(this) AliITStrackerSA(source);
208   return *this;
209  
210 }
211
212 //____________________________________________________________________________
213 AliITStrackerSA::~AliITStrackerSA(){
214   // destructor
215   // if fVertexer is not null, the AliESDVertex obj. is owned by this class
216   // and is deleted here
217   if(fVertexer){
218     if(fVert)delete fVert;
219   }
220   fVert = 0;
221   fVertexer = 0;
222  
223   if(fPhiWin)delete []fPhiWin;
224   if(fLambdaWin)delete []fLambdaWin;
225   fListOfTracks->Delete();
226   delete fListOfTracks;
227   fListOfSATracks->Delete();
228   delete fListOfSATracks;
229   if(fCluLayer){
230     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
231       if(fCluLayer[i]){
232         fCluLayer[i]->Delete();
233         delete fCluLayer[i];
234       }
235     }
236     delete [] fCluLayer;
237   }
238   if(fCluCoord){
239     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
240       if(fCluCoord[i]){
241         fCluCoord[i]->Delete();
242         delete fCluCoord[i];
243       }
244     }
245     delete [] fCluCoord;
246   }
247   
248 }
249
250 //____________________________________________________________________________
251 Int_t AliITStrackerSA::Clusters2Tracks(AliESDEvent *event){
252 // This method is used to find and fit the tracks. By default the corresponding
253 // method in the parent class is invoked. In this way a combined tracking
254 // TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
255 // is done in the ITS only. In the standard reconstruction chain this option
256 // can be set via AliReconstruction::SetOption("ITS","onlyITS")
257   Int_t rc=0;
258   if(!fITSStandAlone){
259     rc=AliITStrackerMI::Clusters2Tracks(event);
260   }
261   else {
262     AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n");
263   }
264   if(!rc) rc=FindTracks(event);
265   return rc;
266 }
267
268 //____________________________________________________________________________
269 void AliITStrackerSA::Init(){
270   //  Reset all data members
271     fPhiEstimate=0;
272     for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
273     fLambdac=0;
274     fPhic=0;
275     fCoef1=0;
276     fCoef2=0;
277     fCoef3=0;
278     fPointc[0]=0;
279     fPointc[1]=0;
280     fVert = 0;
281     fVertexer = 0;
282     Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA();
283     if(nLoops==33){
284       SetFixedWindowSizes();
285     }else{
286       Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA();
287       Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA();
288       Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA();
289       Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA();
290       SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax);
291     }
292     fITSclusters = 0;
293     SetSixPoints();
294     SetOuterStartLayer(0);
295     SetSAFlag(kFALSE);
296     fListOfTracks=new TObjArray(0,0);
297     fListOfSATracks=new TObjArray(0,0);
298     fCluLayer = 0;
299     fCluCoord = 0;
300  }
301 //_______________________________________________________________________
302 void AliITStrackerSA::ResetForFinding(){
303   //  Reset data members used in all loops during track finding
304     fPhiEstimate=0;
305     for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
306     fLambdac=0;
307     fPhic=0;
308     fCoef1=0;
309     fCoef2=0;
310     fCoef3=0;
311     fPointc[0]=0;
312     fPointc[1]=0;
313     fListOfTracks->Delete();
314     fListOfSATracks->Delete();
315 }
316
317  
318
319 //______________________________________________________________________
320 Int_t AliITStrackerSA::FindTracks(AliESDEvent* event){
321
322 // Track finder using the ESD object
323
324
325   
326   if(!fITSclusters){
327     Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
328     return -1;
329   }
330   //Reads event and mark clusters of traks already found, with flag kITSin
331   Int_t nentr=event->GetNumberOfTracks();
332   if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kFALSE) {
333     while (nentr--) {
334       AliESDtrack *track=event->GetTrack(nentr);
335       if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
336         Int_t idx[12];
337         Int_t ncl = track->GetITSclusters(idx);
338         for(Int_t k=0;k<ncl;k++){
339           AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
340           cll->SetBit(kSAflag);
341         }
342       }
343     }
344   }
345   //Get primary vertex
346   Double_t primaryVertex[3];
347   event->GetVertex()->GetXYZ(primaryVertex);
348   //Creates TClonesArray with clusters for each layer. The clusters already used
349   //by AliITStrackerMI are not considered
350   Int_t nclusters[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
351   Int_t dmar[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
352   if (fCluLayer == 0) {
353     fCluLayer = new TClonesArray*[AliITSgeomTGeo::kNLayers];
354     fCluCoord = new TClonesArray*[AliITSgeomTGeo::kNLayers];
355     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
356       fCluLayer[i]=0;
357       fCluCoord[i]=0;
358     }
359   }
360   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
361     AliITSlayer &layer=fgLayers[i];
362     if (!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i)) {
363       for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
364         AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
365         if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
366         if(cls->GetQ()==0) continue; //fake clusters dead zones
367         nclusters[i]++;
368       }
369     }
370     dmar[i]=0;
371     delete fCluLayer[i];
372     fCluLayer[i] = new TClonesArray("AliITSRecPoint",nclusters[i]);
373     delete fCluCoord[i];
374     fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
375   }
376
377   for(Int_t ilay=0;ilay<AliITSgeomTGeo::GetNLayers();ilay++){
378     TClonesArray &clulay = *fCluLayer[ilay];
379     TClonesArray &clucoo = *fCluCoord[ilay];
380     AliITSlayer &layer=fgLayers[ilay];
381     if (!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilay)) {
382       for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
383         AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
384         if(cls->TestBit(kSAflag)==kTRUE) continue;
385         if(cls->GetQ()==0) continue;
386         Double_t phi=0;Double_t lambda=0;
387         Float_t x=0;Float_t y=0;Float_t z=0;
388         Float_t sx=0;Float_t sy=0;Float_t sz=0;
389         GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
390         GetCoorErrors(cls,sx,sy,sz);
391         new (clulay[dmar[ilay]]) AliITSRecPoint(*cls);
392         new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
393         dmar[ilay]++;
394       }
395     }
396   }
397    
398   // track counter
399   Int_t ntrack=0;
400
401   Int_t minNPoints = (fSixPoints ? AliITSgeomTGeo::GetNLayers() : AliITSgeomTGeo::GetNLayers()-1);  
402   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
403     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i)) {
404       minNPoints--;
405     }
406   }
407
408   static Int_t nClusLay[AliITSgeomTGeo::kNLayers];//counter for clusters on each layer
409
410   //loop on different minNPoints
411   Int_t minMinNPoints=minNPoints;
412   if(AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad()) minMinNPoints=2;
413   for(Int_t iMinNPoints=minNPoints; iMinNPoints>=minMinNPoints; iMinNPoints--) {
414     //loop on the different windows
415     for(Int_t nloop=0;nloop<fNloop;nloop++){
416       for(Int_t ncl=0;ncl<fCluLayer[0]->GetEntries();ncl++){ //loop starting from layer 0
417       
418         ResetForFinding();
419         Int_t pflag=0;
420       
421         AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[0]->At(ncl);
422       
423         if(!cl) continue;
424         if (cl->GetQ()<=0) continue;
425       
426         AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(0,ncl); 
427         fPhic = arr->GetPhi();
428         fLambdac = arr->GetLambda();
429         if (TMath::Abs(fLambdac)>0.26*TMath::Pi()) continue;
430         fPhiEstimate = fPhic;
431         AliITStrackSA* trs = new AliITStrackSA(); 
432         fPoint1[0]=primaryVertex[0];
433         fPoint1[1]=primaryVertex[1];
434         
435         
436         fPoint2[0]=arr->GetX();
437         fPoint2[1]=arr->GetY();
438         for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) nClusLay[i]=0;
439         nClusLay[0] = SearchClusters(0,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
440         nClusLay[1] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
441         if(nClusLay[1]>0){
442           pflag=1;
443           fPoint3[0] = fPointc[0];
444           fPoint3[1] = fPointc[1];
445         }
446         nClusLay[2] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
447         if(nClusLay[1]==0 && nClusLay[2]==0) pflag=0;
448         if(nClusLay[2]!=0 && nClusLay[1]!=0){ pflag=1; UpdatePoints();}
449         if(nClusLay[2]!=0 && nClusLay[1]==0){
450           pflag=1;
451           fPoint3[0]=fPointc[0];
452           fPoint3[1]=fPointc[1];
453         }
454         
455         nClusLay[3] = SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
456         pflag=1;
457         if(nClusLay[3]!=0) UpdatePoints();
458         nClusLay[4] = SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag); 
459         pflag=1;
460         if(nClusLay[4]!=0) UpdatePoints();
461         nClusLay[5] = SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag); 
462         
463         
464         Int_t layOK=0;
465         //check of the candidate track
466         for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers();nnp++) {
467           if(nClusLay[nnp]!=0) layOK+=1;
468         }
469         
470         if(layOK>=iMinNPoints) {
471           //printf("-NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]);
472           AliITStrackV2* tr2 = 0;
473           tr2 = FitTrack(trs,primaryVertex);
474           if(!tr2) continue;
475           //printf("-NPOINTS fit: %d\n",tr2->GetNumberOfClusters());      
476
477           
478           AliESDtrack outtrack;
479           outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
480           event->AddTrack(&outtrack);
481           ntrack++;
482         }
483         delete trs;
484       }//end loop on clusters of layer1
485     }//end loop on windows
486   }//end loop on min points
487
488
489   minNPoints--;
490
491   //if 5/6 points are required, second loop starting 
492   //from second layer (SPD2), to find tracks with point of 
493   //layer 1 missing. Not done for cosmics.   
494   if(!fSixPoints && fOuterStartLayer==0) {
495     //printf("looking from SPD2\n");
496     //   counter for clusters on each layer  
497     for(Int_t nloop=0;nloop<fNloop;nloop++){
498       Int_t ncl2=fCluLayer[1]->GetEntries();
499       while(ncl2--){ //loop starting from layer 2
500         ResetForFinding();
501         Int_t pflag=0;
502         AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[1]->At(ncl2);
503         
504         if(!cl) continue;
505         AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(1,ncl2);
506         fPhic = arr->GetPhi();
507         fLambdac = arr->GetLambda();
508         fPhiEstimate = fPhic;
509         
510         AliITStrackSA* trs = new AliITStrackSA(); 
511         fPoint1[0]=primaryVertex[0];
512         fPoint1[1]=primaryVertex[1];
513         
514         fPoint2[0]=arr->GetX();
515         fPoint2[1]=arr->GetY();
516         for(Int_t kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
517         nClusLay[0] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],
518                                trs,primaryVertex[2],pflag);
519         nClusLay[1] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],
520                                trs,primaryVertex[2],pflag);
521         if(nClusLay[1]!=0){
522           pflag=1;
523           fPoint3[0]=fPointc[0];
524           fPoint3[1]=fPointc[1];
525         }
526         nClusLay[2]= SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],
527                               trs,primaryVertex[2],pflag);
528         if(nClusLay[2]!=0){
529           pflag=1;
530           UpdatePoints();
531         }
532         nClusLay[3]= SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],
533                               trs,primaryVertex[2],pflag);
534         if(nClusLay[3]!=0){
535           pflag=1;
536           UpdatePoints();
537         }
538         nClusLay[4]=SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],
539                              trs,primaryVertex[2],pflag);
540         
541         Int_t layOK=0;
542         for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-1;nnp++){
543           if(nClusLay[nnp]!=0) layOK+=1;
544         }
545         if(layOK>=minNPoints){  // 5/6
546           //printf("--NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]);
547        
548           AliITStrackV2* tr2 = 0;
549           tr2 = FitTrack(trs,primaryVertex);
550           if(!tr2) continue;
551           //printf("--NPOINTS fit: %d\n",tr2->GetNumberOfClusters());
552           
553           AliESDtrack outtrack;
554           outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
555           event->AddTrack(&outtrack);
556           ntrack++;
557           
558         }   
559         
560         delete trs;
561       }//end loop on clusters of layer2
562     }
563     
564   } //end opt="5/6"  
565
566   // search for tracks starting from SPD2, SDD1, SDD2, SSD2
567   // for cosmics (A. Dainese 31.07.07)
568   if(fOuterStartLayer>0 && !AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad()) {
569     for(Int_t innLay=1; innLay<=fOuterStartLayer; innLay++) {
570       //printf("Searching from layer %d outward\n",innLay);
571       minNPoints=AliITSgeomTGeo::GetNLayers()-innLay;
572       for(Int_t i=innLay;i<AliITSgeomTGeo::GetNLayers();i++) 
573         if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i)) 
574           minNPoints--;
575       //   counter for clusters on each layer  
576
577       for(Int_t nloop=0;nloop<fNloop;nloop++){
578         Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
579         while(nclInnLay--){ //loop starting from layer innLay
580           ResetForFinding();
581           Int_t pflag=0;
582           AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[innLay]->At(nclInnLay);
583           
584           if(!cl) continue;
585           AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(innLay,nclInnLay);
586           fPhic = arr->GetPhi();
587           fLambdac = arr->GetLambda();
588           fPhiEstimate = fPhic;
589           
590           AliITStrackSA* trs = new AliITStrackSA(); 
591           fPoint1[0]=primaryVertex[0];
592           fPoint1[1]=primaryVertex[1];
593           fPoint2[0]=arr->GetX();
594           fPoint2[1]=arr->GetY();
595           
596           Int_t kk;
597           for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
598           
599           kk=0;
600           nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
601                                   trs,primaryVertex[2],pflag);
602           for(Int_t nextLay=innLay+1; nextLay<AliITSgeomTGeo::GetNLayers(); nextLay++) {
603             kk++;
604             nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
605                                     trs,primaryVertex[2],pflag);
606             if(nClusLay[kk]!=0){
607               pflag=1;
608               if(kk==1) {
609                 fPoint3[0]=fPointc[0];
610                 fPoint3[1]=fPointc[1];
611               } else {
612                 UpdatePoints();
613               }
614             }
615           }
616           
617           Int_t layOK=0;
618           for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
619             if(nClusLay[nnp]!=0) layOK+=1;
620           }
621           if(layOK>=minNPoints){ 
622             //printf("---NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]);
623             AliITStrackV2* tr2 = 0;
624             tr2 = FitTrack(trs,primaryVertex);
625             if(!tr2) continue;
626             //printf("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters());
627             
628             AliESDtrack outtrack;
629             outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
630             event->AddTrack(&outtrack);
631             ntrack++;
632             
633           }   
634           
635           delete trs;
636         }//end loop on clusters of innLay
637       } //end loop on window sizes
638       
639     } //end loop on innLay
640   } //end if(fOuterStartLayer>0)
641   
642
643   // search for 1-point tracks in SPD, only for cosmics
644   // (A.Dainese 21.03.08)
645   if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() && 
646      TMath::Abs(event->GetMagneticField())<0.01) {
647     Int_t outerLayer=1; // only SPD
648     for(Int_t innLay=0; innLay<=TMath::Min(1,fOuterStartLayer); innLay++) {
649       //   counter for clusters on each layer  
650
651       for(Int_t nloop=0;nloop<fNloop;nloop++){
652         Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
653         while(nclInnLay--){ //loop starting from layer innLay
654           ResetForFinding();
655           Int_t pflag=0;
656           AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[innLay]->At(nclInnLay);
657           
658           if(!cl) continue;
659           AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(innLay,nclInnLay);
660           fPhic = arr->GetPhi();
661           fLambdac = arr->GetLambda();
662           fPhiEstimate = fPhic;
663           
664           AliITStrackSA* trs = new AliITStrackSA(); 
665           fPoint1[0]=primaryVertex[0];
666           fPoint1[1]=primaryVertex[1];
667           fPoint2[0]=arr->GetX();
668           fPoint2[1]=arr->GetY();
669           
670           Int_t kk;
671           for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
672           
673           kk=0;
674           nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
675                                   trs,primaryVertex[2],pflag);
676           for(Int_t nextLay=innLay+1; nextLay<=outerLayer; nextLay++) {
677             kk++;
678             nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
679                                     trs,primaryVertex[2],pflag);
680             if(nClusLay[kk]!=0){
681               pflag=1;
682               if(kk==1) {
683                 fPoint3[0]=fPointc[0];
684                 fPoint3[1]=fPointc[1];
685               } else {
686                 UpdatePoints();
687               }
688             }
689           }
690           
691           Int_t layOK=0;
692           for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
693             if(nClusLay[nnp]!=0) layOK+=1;
694           }
695           if(layOK==1) {
696             //printf("----NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]);
697             AliITStrackV2* tr2 = 0;
698             Bool_t onePoint = kTRUE;
699             tr2 = FitTrack(trs,primaryVertex,onePoint);
700             if(!tr2) continue;
701             //printf("----NPOINTS fit: %d\n",tr2->GetNumberOfClusters());
702             
703             AliESDtrack outtrack;
704             outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
705             event->AddTrack(&outtrack);
706             ntrack++;
707             
708           }   
709           
710           delete trs;
711         }//end loop on clusters of innLay
712       } //end loop on window sizes
713       
714     } //end loop on innLay
715   } // end search 1-point tracks
716   
717   Info("FindTracks","Number of found tracks: %d",event->GetNumberOfTracks());
718   return 0;
719
720 }
721  
722 //________________________________________________________________________
723
724 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) {
725   //fit of the found track (most general case, also <6 points, layers missing)
726   // A.Dainese 16.11.07 
727
728   
729   const Int_t kMaxClu=AliITStrackSA::kMaxNumberOfClusters;
730
731   static Int_t firstmod[AliITSgeomTGeo::kNLayers];
732   
733   static Int_t clind[AliITSgeomTGeo::kNLayers][kMaxClu];
734   static Int_t clmark[AliITSgeomTGeo::kNLayers][kMaxClu];
735   static Int_t nnn[AliITSgeomTGeo::kNLayers];
736   static Int_t kkk[AliITSgeomTGeo::kNLayers];
737   static Int_t end[AliITSgeomTGeo::kNLayers];
738   static AliITSRecPoint *listlayer[AliITSgeomTGeo::kNLayers][kMaxClu];
739
740   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
741     firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
742     nnn[i]=0;
743     kkk[i]=0;
744     end[i]=0;
745     for(Int_t j=0;j<kMaxClu; j++){
746       clind[i][j]=0;
747       clmark[i][j]=0;
748       listlayer[i][j]=0;
749    }
750   }
751   
752
753   Int_t nclusters = tr->GetNumberOfClustersSA();
754   for(Int_t ncl=0;ncl<nclusters;ncl++){
755     Int_t index = tr->GetClusterIndexSA(ncl); 
756     AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
757     if(cl->TestBit(kSAflag)==kTRUE) cl->ResetBit(kSAflag);
758     Int_t lay = (index & 0xf0000000) >> 28;
759     Int_t nInLay=end[lay];
760     listlayer[lay][nInLay]=cl;
761     Int_t ind=nnn[lay];
762     clind[lay][ind]=index;
763     end[lay]++;
764   }
765
766   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
767     for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
768       Int_t mark = tr->GetClusterMark(nlay,ncl);
769       Int_t ind=kkk[nlay];
770       clmark[nlay][ind]=mark;
771     }
772   }
773
774
775   Int_t firstLay=-1,secondLay=-1;
776   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
777     if(end[i]==0) {
778       end[i]=1;
779     }else{
780       if(firstLay==-1) {
781         firstLay=i;
782       } else if(secondLay==-1) {
783         secondLay=i;
784       }
785     }
786   }
787
788   if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
789   
790   for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
791     AliITSRecPoint* cl0 = (AliITSRecPoint*)listlayer[0][l0];
792     for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
793       AliITSRecPoint* cl1 = (AliITSRecPoint*)listlayer[1][l1];
794       for(Int_t l2=0;l2<end[2];l2++){  //loop on layer 3
795         AliITSRecPoint* cl2 = (AliITSRecPoint*)listlayer[2][l2];
796         for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4   
797           AliITSRecPoint* cl3 = (AliITSRecPoint*)listlayer[3][l3];
798           for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
799             AliITSRecPoint* cl4 = (AliITSRecPoint*)listlayer[4][l4];
800             for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6  
801               AliITSRecPoint* cl5 = (AliITSRecPoint*)listlayer[5][l5];
802
803
804               Double_t x1,y1,z1,sx1,sy1,sz1;
805               Double_t x2,y2,z2,sx2,sy2,sz2;
806               AliITSRecPoint* p1=0;
807               AliITSRecPoint* p2=0;
808               Int_t index1=0,index2=0;
809               Int_t mrk1=0,mrk2=0;
810
811               switch(firstLay) {
812               case 0:
813                 p1=cl0;
814                 index1=clind[0][l0];mrk1=clmark[0][l0];
815                 break;
816               case 1:
817                 p1=cl1;
818                 index1=clind[1][l1];mrk1=clmark[1][l1];
819                 break;
820               case 2:
821                 p1=cl2;
822                 index1=clind[2][l2];mrk1=clmark[2][l2];
823                 break;
824               case 3:
825                 p1=cl3;
826                 index1=clind[3][l3];mrk1=clmark[3][l3];
827                 break;
828               case 4:
829                 p1=cl4;
830                 index1=clind[4][l4];mrk1=clmark[4][l4];
831                 break;
832               }
833
834               switch(secondLay) {
835               case 1:
836                 p2=cl1;
837                 index2=clind[1][l1];mrk2=clmark[1][l1];
838                 break;
839               case 2:
840                 p2=cl2;
841                 index2=clind[2][l2];mrk2=clmark[2][l2];
842                 break;
843               case 3:
844                 p2=cl3;
845                 index2=clind[3][l3];mrk2=clmark[3][l3];
846                 break;
847               case 4:
848                 p2=cl4;
849                 index2=clind[4][l4];mrk2=clmark[4][l4];
850                 break;
851               case 5:
852                 p2=cl5;
853                 index2=clind[5][l5];mrk2=clmark[5][l5];
854                 break;
855               default:
856                 p2=0;
857                 index2=-1;mrk2=-1;
858                 break;
859               }
860
861               Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay]; 
862               Int_t layer,ladder,detector;
863               AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
864               Float_t yclu1 = p1->GetY();
865               Float_t zclu1 = p1->GetZ();
866               Double_t cv=0,tgl2=0,phi2=0;
867               
868               Int_t cln1=mrk1;
869               AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,cln1);
870               x1 = arr1->GetX();
871               y1 = arr1->GetY();
872               z1 = arr1->GetZ();
873               sx1 = arr1->GetSx();
874               sy1 = arr1->GetSy();
875               sz1 = arr1->GetSz();
876
877               if(secondLay>0) {
878                 Int_t cln2=mrk2;
879                 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,cln2);
880                 x2 = arr2->GetX();
881                 y2 = arr2->GetY();
882                 z2 = arr2->GetZ();
883                 sx2 = arr2->GetSx();
884                 sy2 = arr2->GetSy();
885                 sz2 = arr2->GetSz();
886                 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
887                 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
888                 phi2 = TMath::ATan2((y2-y1),(x2-x1));
889               } else { // special case of 1-point tracks, only for cosmics (B=0)
890                 x2 = primaryVertex[0];
891                 y2 = primaryVertex[1];
892                 z2 = primaryVertex[2];
893                 cv = 0;
894                 tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
895                 phi2 = TMath::ATan2((y1-y2),(x1-x2));
896               }
897
898
899               AliITStrackSA* trac = new AliITStrackSA(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
900
901
902               if(cl5!=0) {
903                 trac->AddClusterV2(5,(clind[5][l5] & 0x0fffffff)>>0);
904                 trac->AddClusterMark(5,clmark[5][l5]);
905               }
906               if(cl4!=0){
907                 trac->AddClusterV2(4,(clind[4][l4] & 0x0fffffff)>>0);
908                 trac->AddClusterMark(4,clmark[4][l4]);
909               }
910               if(cl3!=0){
911                 trac->AddClusterV2(3,(clind[3][l3] & 0x0fffffff)>>0);
912                 trac->AddClusterMark(3,clmark[3][l3]);
913               }
914               if(cl2!=0){
915                 trac->AddClusterV2(2,(clind[2][l2] & 0x0fffffff)>>0);
916                 trac->AddClusterMark(2,clmark[2][l2]);
917               }
918               if(cl1!=0){
919                 trac->AddClusterV2(1,(clind[1][l1] & 0x0fffffff)>>0);
920                 trac->AddClusterMark(1,clmark[1][l1]);
921               }
922               if(cl0!=0){
923                 trac->AddClusterV2(0,(clind[0][l0] & 0x0fffffff)>>0);
924                 trac->AddClusterMark(0,clmark[0][l0]);
925               }
926
927               //fit with Kalman filter using AliITStrackerMI::RefitAt()
928               AliITStrackMI* ot = new AliITStrackSA(*trac);
929               
930               ot->ResetCovariance(10.);
931               ot->ResetClusters();
932               
933               if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),ot,trac)){ //fit from layer 1 to layer 6
934                 AliITStrackMI *otrack2 = new AliITStrackMI(*ot);
935                 otrack2->ResetCovariance(10.); 
936                 otrack2->ResetClusters();
937                 //fit from layer 6 to layer 1
938                 if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),otrack2,ot)) {
939                   fListOfTracks->AddLast(otrack2);
940                   fListOfSATracks->AddLast(trac);
941                 } else {
942                   delete otrack2;
943                   delete trac;
944                 }
945                               
946               }       
947           
948               delete ot;
949             }//end loop layer 6
950           }//end loop layer 5
951         }//end loop layer 4        
952       }//end loop layer 3
953     }//end loop layer 2 
954   }//end loop layer 1
955
956
957
958
959   if(fListOfTracks->GetEntries()==0) return 0;
960
961   Int_t lowchi2 = FindTrackLowChiSquare();
962   AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
963   AliITStrackSA* trsa = (AliITStrackSA*)fListOfSATracks->At(lowchi2);
964  
965   if(otrack==0) return 0;
966
967   Int_t indexc[AliITSgeomTGeo::kNLayers];
968   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) indexc[i]=0;
969   for(Int_t nind=0;nind<otrack->GetNumberOfClusters();nind++){
970     indexc[nind] = otrack->GetClusterIndex(nind);
971   }      
972   Int_t labl[6][3];
973   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
974     if(i<otrack->GetNumberOfClusters()) {
975       AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc[i]);
976       labl[i][0]=cl->GetLabel(0);
977       labl[i][1]=cl->GetLabel(1);
978       labl[i][2]=cl->GetLabel(2);
979     } else {
980       labl[i][0]=-1;
981       labl[i][1]=-1;
982       labl[i][2]=-1;
983     }
984   }
985
986   CookLabel(otrack,0.); //MI change - to see fake ratio
987  
988   Int_t label=FindLabel(labl[0][0],labl[1][0],labl[2][0],labl[3][0],labl[4][0],labl[5][0]);
989   Int_t lflag=0;
990   for(Int_t i=0;i<otrack->GetNumberOfClusters();i++)
991     if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
992   
993   if(lflag<otrack->GetNumberOfClusters()) label = -label;
994   otrack->SetLabel(label);  
995
996   //remove clusters of found track
997   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
998     for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
999       Int_t index = trsa->GetClusterMark(nlay,cln);
1000       fCluLayer[nlay]->RemoveAt(index);
1001       RemoveClusterCoord(nlay,index);
1002       fCluLayer[nlay]->Compress();
1003     }    
1004   }
1005
1006   return otrack;
1007
1008 }
1009
1010
1011
1012 //_______________________________________________________
1013 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
1014   //function used to to find the clusters associated to the track
1015
1016   if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(layer)) return 0;
1017
1018   Int_t nc=0;
1019   AliITSlayer &lay = fgLayers[layer];
1020   Double_t r=lay.GetR();
1021   if(pflag==1){      
1022     Float_t cx1,cx2,cy1,cy2;
1023     FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
1024     if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
1025        return 0;
1026     Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
1027     Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
1028     fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
1029   }
1030
1031  
1032   Int_t ncl = fCluLayer[layer]->GetEntries();
1033   for (Int_t index=0; index<ncl; index++) {
1034     AliITSRecPoint *c = (AliITSRecPoint*)fCluLayer[layer]->At(index);
1035     if (!c) continue;
1036     if (c->GetQ()<=0) continue;
1037     
1038      AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
1039      Double_t phi = arr->GetPhi();
1040      if (TMath::Abs(phi-fPhiEstimate)>phiwindow) continue;
1041
1042      Double_t lambda = arr->GetLambda();
1043      if (TMath::Abs(lambda-fLambdac)>lambdawindow) continue;
1044
1045      if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
1046      if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
1047      Int_t orind = arr->GetOrInd();
1048      trs->AddClusterSA(layer,orind);
1049      trs->AddClusterMark(layer,index);
1050        
1051      nc++;
1052      fLambdac=lambda;
1053      fPhiEstimate=phi;
1054
1055      fPointc[0]=arr->GetX();
1056      fPointc[1]=arr->GetY();
1057
1058   }
1059   return nc;
1060 }
1061
1062 //________________________________________________________________
1063 void AliITStrackerSA::UpdatePoints(){
1064   //update of points for the estimation of the curvature  
1065
1066   fPoint2[0]=fPoint3[0];
1067   fPoint2[1]=fPoint3[1];
1068   fPoint3[0]=fPointc[0];
1069   fPoint3[1]=fPointc[1];
1070
1071   
1072 }
1073
1074 //___________________________________________________________________
1075 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){
1076
1077    //given (x,y) of three recpoints (in global coordinates) 
1078    //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
1079
1080    Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
1081    if(den==0) return 0;
1082    a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
1083    b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
1084    c = -x1*x1-y1*y1-a*x1-b*y1;
1085    return 1;
1086  }
1087 //__________________________________________________________________________
1088  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){
1089  
1090  //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
1091  //c2 is -rlayer*rlayer
1092
1093   if(a1==0) return 0;
1094  Double_t m = c2-c1; 
1095  Double_t aA = (b1*b1)/(a1*a1)+1;
1096  Double_t bB = (-2*m*b1/(a1*a1));
1097  Double_t cC = c2+(m*m)/(a1*a1);
1098  Double_t dD = bB*bB-4*aA*cC;
1099  if(dD<0) return 0;
1100  
1101  y1 = (-bB+TMath::Sqrt(dD))/(2*aA); 
1102  y2 = (-bB-TMath::Sqrt(dD))/(2*aA); 
1103  x1 = (c2-c1-b1*y1)/a1;
1104  x2 = (c2-c1-b1*y2)/a1;
1105
1106  return 1; 
1107 }
1108 //____________________________________________________________________
1109 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t 
1110 x2,Double_t y2,Double_t x3,Double_t y3){
1111
1112   //calculates the curvature of track  
1113   Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
1114   if(den==0) return 0;
1115   Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
1116   Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
1117   Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
1118   Double_t xc=-a/2.;
1119
1120   if((a*a+b*b-4*c)<0) return 0;
1121   Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
1122   if(rad==0) return 0;
1123   
1124   if((x1>0 && y1>0 && x1<xc)) rad*=-1;
1125   if((x1<0 && y1>0 && x1<xc)) rad*=-1;
1126   //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
1127   // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
1128   
1129   return 1/rad;
1130  
1131 }
1132
1133
1134 //____________________________________________________________________
1135 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
1136
1137   //Returns the point closest to pp
1138
1139   Double_t diff1 = p1-pp;
1140   Double_t diff2 = p2-pp;
1141   
1142   if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
1143   else fPhiEstimate=p2;  
1144   return fPhiEstimate;
1145   
1146 }
1147
1148
1149 //_________________________________________________________________
1150 Int_t AliITStrackerSA::FindTrackLowChiSquare() const {
1151   // returns track with lowest chi square  
1152   Int_t dim=fListOfTracks->GetEntries();
1153   if(dim<=1) return 0;
1154   AliITStrackV2* trk = (AliITStrackV2*)fListOfTracks->At(0);
1155   Double_t minChi2=trk->GetChi2();
1156   Int_t index=0;
1157   for(Int_t i=1;i<dim;i++){
1158     trk = (AliITStrackV2*)fListOfTracks->At(i);
1159     Double_t chi2=trk->GetChi2();
1160     if(chi2<minChi2){
1161       minChi2=chi2;
1162       index=i;
1163     }
1164   }
1165   return index;
1166 }
1167
1168 //__________________________________________________________
1169 Int_t AliITStrackerSA::FindLabel(Int_t l0, Int_t l1, Int_t l2, Int_t l3, Int_t l4, Int_t l5){
1170
1171   //function used to determine the track label
1172   
1173   Int_t lb[6] = {l0,l1,l2,l3,l4,l5};
1174   Int_t aa[6]={1,1,1,1,1,1};
1175   Int_t ff=0; 
1176   Int_t ll=0;
1177   Int_t k=0;Int_t w=0;Int_t num=6;
1178   for(Int_t i=5;i>=0;i--) if(lb[i]==-1) num=i;
1179   
1180   while(k<num){
1181   
1182     for(Int_t i=k+1;i<num;i++){
1183     
1184       if(lb[k]==lb[i] && aa[k]!=0){
1185       
1186         aa[k]+=1;
1187         aa[i]=0;
1188       }
1189     }
1190     k++;
1191   }
1192
1193   while(w<num){
1194  
1195     for(Int_t j=0;j<6;j++){
1196       if(aa[w]<aa[j]){
1197       ff=aa[w];
1198       aa[w]=aa[j];
1199       aa[j]=ff;
1200       ll=lb[w];
1201       lb[w]=lb[j];
1202       lb[j]=ll;
1203      }
1204     }
1205     w++;
1206   }
1207   
1208   if(num<1) return -1; 
1209   return lb[num-1];
1210 }
1211
1212 //_____________________________________________________________________________
1213 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,
1214 Int_t gl12, Int_t gl13, Int_t gl14,Int_t gl15, Int_t gl16, Int_t gl17, Int_t gl18, Int_t minNPoints){
1215
1216  
1217   //function used to assign label to the found track. If track is fake, the label is negative
1218
1219   Int_t lb0[6] = {gl1,gl2,gl3,gl4,gl5,gl6};
1220   Int_t lb1[6] = {gl7,gl8,gl9,gl10,gl11,gl12};
1221   Int_t lb2[6] = {gl13,gl14,gl15,gl16,gl17,gl18};
1222   Int_t ll=FindLabel(lb0[0],lb0[1],lb0[2],lb0[3],lb0[4],lb0[5]);
1223   Int_t lflag=0;Int_t num=6;
1224   if(lb0[5]==-1 && lb1[5]==-1 && lb2[5]==-1) num=5;
1225
1226   for(Int_t i=0;i<num;i++){
1227     if(lb0[i]==ll || lb1[i]==ll || lb2[i]==ll) lflag+=1;
1228   }
1229
1230   if(lflag>=minNPoints) return ll;
1231   else return -ll;
1232
1233   
1234 }
1235 //_____________________________________________________________________________
1236 void AliITStrackerSA::SetCalculatedWindowSizes(Int_t n, Float_t phimin, Float_t phimax, Float_t lambdamin, Float_t lambdamax){
1237   // Set sizes of the phi and lambda windows used for track finding
1238   fNloop = n;
1239   if(fPhiWin) delete [] fPhiWin;
1240   if(fLambdaWin) delete [] fLambdaWin;
1241   fPhiWin = new Double_t[fNloop];
1242   fLambdaWin = new Double_t[fNloop];
1243   Float_t stepPhi=(phimax-phimin)/(Float_t)(fNloop-1);
1244   Float_t stepLambda=(lambdamax-lambdamin)/(Float_t)(fNloop-1);
1245   for(Int_t k=0;k<fNloop;k++){
1246     Float_t phi=phimin+k*stepPhi;
1247     Float_t lam=lambdamin+k*stepLambda;
1248     fPhiWin[k]=phi;
1249     fLambdaWin[k]=lam;
1250   }
1251 }
1252 //_____________________________________________________________________________
1253 void AliITStrackerSA::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
1254   // Set sizes of the phi and lambda windows used for track finding
1255   fNloop = n;
1256   if(phi){ // user defined values
1257     fPhiWin = new Double_t[fNloop];
1258     fLambdaWin = new Double_t[fNloop];
1259     for(Int_t k=0;k<fNloop;k++){
1260       fPhiWin[k]=phi[k];
1261       fLambdaWin[k]=lam[k];
1262     }
1263   }
1264   else {  // default values
1265             
1266     Double_t phid[33]   = {0.002,0.003,0.004,0.0045,0.0047,
1267                            0.005,0.0053,0.0055,
1268                            0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
1269                            0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
1270                            0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
1271     Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
1272                             0.005,0.005,0.006,
1273                             0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1274                             0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1275                             0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1276     
1277     if(fNloop!=33){
1278       fNloop = 33;
1279     }
1280     
1281     
1282     fPhiWin = new Double_t[fNloop];
1283     fLambdaWin = new Double_t[fNloop];
1284
1285     Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
1286   
1287     for(Int_t k=0;k<fNloop;k++){
1288       fPhiWin[k]=phid[k]*factor;
1289       fLambdaWin[k]=lambdad[k]*factor;
1290     }
1291   
1292   }
1293
1294 }
1295 //_______________________________________________________________________
1296 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Float_t &x, Float_t &y,Float_t &z,Double_t* vertex){
1297   //Returns values of phi (azimuthal) and lambda angles for a given cluster
1298 /*  
1299   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1300   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1301   Float_t tx,ty,tz;  fGeom->GetTrans(lay,lad,det,tx,ty,tz);     
1302
1303   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1304   Double_t phi1=TMath::Pi()/2+alpha;
1305   if (lay==1) phi1+=TMath::Pi();
1306
1307   Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1308   Float_t r=tx*cp+ty*sp;
1309
1310   xyz= r*cp - cl->GetY()*sp;
1311   y= r*sp + cl->GetY()*cp;
1312   z=cl->GetZ();
1313 */
1314   Float_t xyz[3];
1315   cl->GetGlobalXYZ(xyz);
1316   x=xyz[0];
1317   y=xyz[1];
1318   z=xyz[2];
1319  
1320   phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1321   lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1322 }
1323
1324 //________________________________________________________________________
1325 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Float_t &sx,Float_t &sy, Float_t &sz){
1326
1327   //returns sigmax, y, z of cluster in global coordinates
1328 /*
1329   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1330   Int_t lay,lad,det; 
1331   AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1332  
1333   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1334   Double_t phi=TMath::Pi()/2+alpha;
1335   if (lay==1) phi+=TMath::Pi();
1336
1337   Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1338 */
1339   Float_t covm[6];
1340   cl->GetGlobalCov(covm);
1341   sx=TMath::Sqrt(covm[0]);
1342   sy=TMath::Sqrt(covm[3]);
1343   sz=TMath::Sqrt(covm[5]);
1344 /*
1345   sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1346   sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1347   sz = TMath::Sqrt(cl->GetSigmaZ2());
1348 */
1349 }