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