Added option for different binning of DCAxy axis in THnSparse. Same width for all...
[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 ITS tracker class                        //
20 //  Origin:  Elisabetta Crescio - crescio@to.infn.it     //
21 //  Updated: Francesco Prino    - prino@to.infn.it       //
22 ///////////////////////////////////////////////////////////
23
24 #include <stdlib.h>
25
26 #include <TArrayI.h>
27 #include <TBranch.h>
28 #include <TObjArray.h>
29 #include <TTree.h>
30
31 #include "AliESDEvent.h"
32 #include "AliESDVertex.h"
33 #include "AliESDtrack.h"
34 #include "AliITSVertexer.h"
35 #include "AliITSclusterTable.h"
36 #include "AliITSRecPoint.h"
37 #include "AliITSgeomTGeo.h"
38 #include "AliITStrackSA.h"
39 #include "AliITStrackerSA.h"
40 #include "AliITSReconstructor.h"
41 #include "AliLog.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 fListOfTracks(0),
59 fListOfSATracks(0),
60 fITSclusters(0),
61 fInwardFlag(0),
62 fOuterStartLayer(0),
63 fInnerStartLayer(5),
64 fMinNPoints(0),
65 fMinQ(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 fListOfTracks(0),
84 fListOfSATracks(0),
85 fITSclusters(0),
86 fInwardFlag(0),
87 fOuterStartLayer(0),
88 fInnerStartLayer(5),
89 fMinNPoints(0),
90 fMinQ(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  
100 }
101
102 //____________________________________________________________________________
103 AliITStrackerSA::~AliITStrackerSA(){
104   // destructor
105  
106   if(fPhiWin)delete []fPhiWin;
107   if(fLambdaWin)delete []fLambdaWin;
108   fListOfTracks->Delete();
109   delete fListOfTracks;
110   fListOfSATracks->Delete();
111   delete fListOfSATracks;
112   if(fCluCoord){
113     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
114       if(fCluCoord[i]){
115         fCluCoord[i]->Delete();
116         delete fCluCoord[i];
117       }
118     }
119     delete [] fCluCoord;
120   }
121 }
122
123 //____________________________________________________________________________
124 Int_t AliITStrackerSA::Clusters2Tracks(AliESDEvent *event){
125 // This method is used to find and fit the tracks. By default the corresponding
126 // method in the parent class is invoked. In this way a combined tracking
127 // TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
128 // is done in the ITS only. In the standard reconstruction chain this option
129 // can be set via AliReconstruction::SetOption("ITS","onlyITS")
130   Int_t rc=0;
131
132   if(!fITSStandAlone){
133     rc=AliITStrackerMI::Clusters2Tracks(event);
134   }
135   else {
136     AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n");
137   }
138   if(!rc){ 
139     rc=FindTracks(event,kFALSE);
140     Int_t nSPDcontr=0;
141     const AliESDVertex *spdv = event->GetPrimaryVertexSPD();
142     if(spdv) nSPDcontr = spdv->GetNContributors();
143     if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kTRUE && 
144        nSPDcontr<=AliITSReconstructor::GetRecoParam()->GetMaxSPDcontrForSAToUseAllClusters()) {
145       rc=FindTracks(event,kTRUE);
146     }
147   }
148   return rc;
149 }
150
151 //____________________________________________________________________________
152 void AliITStrackerSA::Init(){
153   //  Reset all data members
154     fPhiEstimate=0;
155     for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
156     fLambdac=0;
157     fPhic=0;
158     fCoef1=0;
159     fCoef2=0;
160     fCoef3=0;
161     fPointc[0]=0;
162     fPointc[1]=0;
163     Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA();
164     if(nLoops==32){
165       SetFixedWindowSizes();
166     }else{
167       Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA();
168       Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA();
169       Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA();
170       Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA();
171       SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax);
172     }
173     fMinQ=AliITSReconstructor::GetRecoParam()->GetSAMinClusterCharge();
174     fITSclusters = 0;
175     SetOuterStartLayer(1);
176     SetSAFlag(kFALSE);
177     fListOfTracks=new TClonesArray("AliITStrackMI",100);
178     fListOfSATracks=new TClonesArray("AliITStrackSA",100);
179     fCluCoord = 0;
180     fMinNPoints = 3;
181  }
182 //_______________________________________________________________________
183 void AliITStrackerSA::ResetForFinding(){
184   //  Reset data members used in all loops during track finding
185     fPhiEstimate=0;
186     for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
187     fLambdac=0;
188     fPhic=0;
189     fCoef1=0;
190     fCoef2=0;
191     fCoef3=0;
192     fPointc[0]=0;
193     fPointc[1]=0;
194     fListOfTracks->Clear();
195     fListOfSATracks->Clear();
196 }
197
198  
199
200 //______________________________________________________________________
201 Int_t AliITStrackerSA::FindTracks(AliESDEvent* event, Bool_t useAllClusters){
202
203 // Track finder using the ESD object
204
205   AliDebug(2,Form(" field is %f",event->GetMagneticField()));
206   AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
207
208   if(!fITSclusters){
209     Fatal("FindTracks","ITS cluster tree is not accessed!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
210     return -1;
211   }
212   //Reads event and mark clusters of traks already found, with flag kITSin
213   Int_t nentr=event->GetNumberOfTracks();
214   if(!useAllClusters) {
215     while (nentr--) {
216       AliESDtrack *track=event->GetTrack(nentr);
217       if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
218         Int_t idx[12];
219         Int_t ncl = track->GetITSclusters(idx);
220         for(Int_t k=0;k<ncl;k++){
221           AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
222           cll->SetBit(kSAflag);
223         }
224       }
225     }
226   }else{
227     while (nentr--) {
228       AliESDtrack *track=event->GetTrack(nentr);
229       if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
230         Int_t idx[12];
231         Int_t ncl = track->GetITSclusters(idx);
232         for(Int_t k=0;k<ncl;k++){
233           AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
234           cll->ResetBit(kSAflag);
235         }
236       }
237     }
238   }
239   //Get primary vertex
240   Double_t primaryVertex[3];
241   event->GetVertex()->GetXYZ(primaryVertex);
242   //Creates TClonesArray with clusters for each layer. The clusters already used
243   //by AliITStrackerMI are not considered
244   Int_t nclusters[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
245   Int_t dmar[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
246   if (fCluCoord == 0) {
247     fCluCoord = new TClonesArray*[AliITSgeomTGeo::kNLayers];
248     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
249       fCluCoord[i]=0;
250     }
251   }
252   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
253     AliITSlayer &layer=fgLayers[i];
254     if (!ForceSkippingOfLayer(i)) {
255       for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
256         AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
257         if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
258         if(cls->GetQ()==0) continue; //fake clusters dead zones
259         if(i>1 && cls->GetQ()<=fMinQ) continue; // cut on SDD and SSD cluster charge
260         nclusters[i]++;
261       }
262     }
263     dmar[i]=0;
264     if(!fCluCoord[i]){
265       fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
266     }else{
267       fCluCoord[i]->Delete();
268       fCluCoord[i]->Expand(nclusters[i]);
269     }
270   }
271
272   for(Int_t ilay=0;ilay<AliITSgeomTGeo::GetNLayers();ilay++){
273     TClonesArray &clucoo = *fCluCoord[ilay];
274     AliITSlayer &layer=fgLayers[ilay];
275     if (!ForceSkippingOfLayer(ilay)) {
276       for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
277         AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
278         if(cls->TestBit(kSAflag)==kTRUE) continue;
279         if(cls->GetQ()==0) continue;
280         if(ilay>1 && cls->GetQ()<=fMinQ) continue; 
281         Double_t phi=0;Double_t lambda=0;
282         Double_t x=0;Double_t y=0;Double_t z=0;
283         Double_t sx=0;Double_t sy=0;Double_t sz=0;
284         GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
285         GetCoorErrors(cls,sx,sy,sz);
286         new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
287         dmar[ilay]++;
288       }
289     }
290     fCluCoord[ilay]->Sort();
291   }
292    
293   // track counter
294   Int_t ntrack=0;
295
296   static Int_t nClusLay[AliITSgeomTGeo::kNLayers];//counter for clusters on each layer
297   Int_t startLayForSeed=0;
298   Int_t lastLayForSeed=fOuterStartLayer;
299   Int_t nSeedSteps=lastLayForSeed-startLayForSeed;
300   Int_t seedStep=1;
301   if(fInwardFlag){
302     startLayForSeed=AliITSgeomTGeo::GetNLayers()-1;
303     lastLayForSeed=fInnerStartLayer;
304     nSeedSteps=startLayForSeed-lastLayForSeed;
305     seedStep=-1;
306   }
307
308   // loop on minimum number of points
309   for(Int_t iMinNPoints=AliITSgeomTGeo::GetNLayers(); iMinNPoints>=fMinNPoints; iMinNPoints--) {
310
311     // loop on starting layer for track finding 
312     for(Int_t iSeedLay=0; iSeedLay<=nSeedSteps; iSeedLay++) {
313       Int_t theLay=startLayForSeed+iSeedLay*seedStep;
314       if(ForceSkippingOfLayer(theLay)) continue;
315       Int_t minNPoints=iMinNPoints-theLay;
316       if(fInwardFlag) minNPoints=iMinNPoints-(AliITSgeomTGeo::GetNLayers()-1-theLay);
317       for(Int_t i=theLay+1;i<AliITSgeomTGeo::GetNLayers();i++)
318         if(ForceSkippingOfLayer(i)) 
319           minNPoints--;
320       if(minNPoints<fMinNPoints) continue;
321
322       // loop on phi and lambda window size
323       for(Int_t nloop=0;nloop<fNloop;nloop++){
324         Int_t nclTheLay=fCluCoord[theLay]->GetEntries();
325         while(nclTheLay--){ 
326           ResetForFinding();
327           Bool_t useRP=SetFirstPoint(theLay,nclTheLay,primaryVertex);
328           if(!useRP) continue;      
329           AliITStrackSA trs;
330             
331           Int_t pflag=0;            
332           Int_t kk;
333           for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
334             
335           kk=0;
336           nClusLay[kk] = SearchClusters(theLay,fPhiWin[nloop],fLambdaWin[nloop],
337                                         &trs,primaryVertex[2],pflag);
338           Int_t nextLay=theLay+seedStep;
339           Bool_t goon=kTRUE;
340           if(nextLay<0 || nextLay == 6) goon = kFALSE;
341           while(goon){
342             kk++;
343             nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
344                                             &trs,primaryVertex[2],pflag);
345             if(nClusLay[kk]!=0){
346               pflag=1;
347               if(kk==1) {
348                 fPoint3[0]=fPointc[0];
349                 fPoint3[1]=fPointc[1];
350               } else {
351                 UpdatePoints();
352               }
353             }
354             nextLay+=seedStep;
355             if(nextLay<0 || nextLay==6) goon=kFALSE;
356           }
357
358             
359           Int_t layOK=0;
360           if(!fInwardFlag){
361             for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-theLay;nnp++){
362               if(nClusLay[nnp]!=0) layOK+=1;
363             }
364           }else{
365             for(Int_t nnp=theLay; nnp>=0; nnp--){
366               if(nClusLay[nnp]!=0) layOK+=1;
367             }
368           }
369           if(layOK>=minNPoints){ 
370             AliDebug(2,Form("---NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]));
371             AliITStrackV2* tr2 = 0;
372             tr2 = FitTrack(&trs,primaryVertex);
373             if(!tr2){ 
374               continue;
375             }
376             AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
377               
378             StoreTrack(tr2,event,useAllClusters);
379             ntrack++;
380               
381           }   
382           
383         }//end loop on clusters of theLay
384       } //end loop on window sizes
385     } //end loop on theLay
386   }//end loop on min points
387
388   // search for 1-point tracks in SPD, only for cosmics
389   // (A.Dainese 21.03.08)
390   if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() && 
391      TMath::Abs(event->GetMagneticField())<0.01) {
392     Int_t outerLayer=1; // only SPD
393     for(Int_t innLay=0; innLay<=TMath::Min(1,fOuterStartLayer); innLay++) {
394       //   counter for clusters on each layer  
395
396       for(Int_t nloop=0;nloop<fNloop;nloop++){
397         Int_t nclInnLay=fCluCoord[innLay]->GetEntries();
398         while(nclInnLay--){ //loop starting from layer innLay
399           ResetForFinding();
400           Bool_t useRP=SetFirstPoint(innLay,nclInnLay,primaryVertex);
401           if(!useRP) continue;
402           AliITStrackSA trs;
403             
404           Int_t pflag=0;            
405           Int_t kk;
406           for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
407           
408           kk=0;
409           nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
410                                   &trs,primaryVertex[2],pflag);
411           for(Int_t nextLay=innLay+1; nextLay<=outerLayer; nextLay++) {
412             kk++;
413             nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
414                                     &trs,primaryVertex[2],pflag);
415             if(nClusLay[kk]!=0){
416               pflag=1;
417               if(kk==1) {
418                 fPoint3[0]=fPointc[0];
419                 fPoint3[1]=fPointc[1];
420               } else {
421                 UpdatePoints();
422               }
423             }
424           }
425           
426           Int_t layOK=0;
427           for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
428             if(nClusLay[nnp]!=0) layOK+=1;
429           }
430           if(layOK==1) {
431             AliDebug(2,Form("----NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]));
432             AliITStrackV2* tr2 = 0;
433             Bool_t onePoint = kTRUE;
434             tr2 = FitTrack(&trs,primaryVertex,onePoint);
435             if(!tr2){
436               continue;
437             }
438             AliDebug(2,Form("----NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
439             
440             StoreTrack(tr2,event,useAllClusters);
441             ntrack++;
442             
443           }   
444           
445         }//end loop on clusters of innLay
446       } //end loop on window sizes
447       
448     } //end loop on innLay
449   } // end search 1-point tracks
450   
451   if(!useAllClusters) AliInfo(Form("Number of found tracks: %d",event->GetNumberOfTracks()));
452   ResetForFinding();
453   return 0;
454
455 }
456  
457 //________________________________________________________________________
458
459 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) {
460   //fit of the found track (most general case, also <6 points, layers missing)
461   // A.Dainese 16.11.07 
462
463   
464   const Int_t kMaxClu=AliITStrackSA::kMaxNumberOfClusters;
465
466   static Int_t firstmod[AliITSgeomTGeo::kNLayers];  
467   static Int_t clind[AliITSgeomTGeo::kNLayers][kMaxClu];
468   static Int_t clmark[AliITSgeomTGeo::kNLayers][kMaxClu];
469   static Int_t end[AliITSgeomTGeo::kNLayers];
470   static Int_t indices[AliITSgeomTGeo::kNLayers];
471
472   static AliITSRecPoint *listlayer[AliITSgeomTGeo::kNLayers][kMaxClu];
473
474   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
475     firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
476     end[i]=0;
477     for(Int_t j=0;j<kMaxClu; j++){
478       clind[i][j]=0;
479       clmark[i][j]=0;
480       listlayer[i][j]=0;
481    }
482   }
483   
484
485   Int_t nclusters = tr->GetNumberOfClustersSA();
486   for(Int_t ncl=0;ncl<nclusters;ncl++){
487     Int_t index = tr->GetClusterIndexSA(ncl); 
488     AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
489     Int_t lay = (index & 0xf0000000) >> 28;
490     Int_t nInLay=end[lay];
491     listlayer[lay][nInLay]=cl;
492     clind[lay][nInLay]=index;
493     end[lay]++;
494   }
495
496   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
497     for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
498       Int_t mark = tr->GetClusterMark(nlay,ncl);
499       clmark[nlay][ncl]=mark;
500     }
501   }
502
503
504   Int_t firstLay=-1,secondLay=-1;
505   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
506     if(end[i]==0) {
507       end[i]=1;
508     }else{
509       if(firstLay==-1) {
510         firstLay=i;
511       } else if(secondLay==-1) {
512         secondLay=i;
513       }
514     }
515   }
516
517   if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
518   TClonesArray &arrMI= *fListOfTracks;
519   TClonesArray &arrSA= *fListOfSATracks;
520   Int_t nFoundTracks=0;
521
522
523   for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
524     indices[0]=l0;
525     for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
526       indices[1]=l1;
527       for(Int_t l2=0;l2<end[2];l2++){  //loop on layer 3
528         indices[2]=l2;
529         for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4   
530           indices[3]=l3;
531           for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
532             indices[4]=l4;
533             for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6  
534               indices[5]=l5;
535
536               // estimate curvature from 2 innermost points (or innermost point + vertex)
537
538               Int_t iFirstLay=indices[firstLay];
539               Int_t mrk1=clmark[firstLay][iFirstLay];
540
541               AliITSRecPoint* p1=(AliITSRecPoint*)listlayer[firstLay][iFirstLay];
542               Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay]; 
543               Int_t layer,ladder,detector;
544               AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
545               Double_t yclu1 = p1->GetY();
546               Double_t zclu1 = p1->GetZ();
547
548               Double_t x1,y1,z1;
549               Double_t x2,y2,z2;
550               Double_t cv=0,tgl2=0,phi2=0;
551               AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,mrk1);
552               x1 = arr1->GetX();
553               y1 = arr1->GetY();
554               z1 = arr1->GetZ();
555
556               if(secondLay>0) {
557                 Int_t iSecondLay=indices[secondLay];          
558                 Int_t mrk2=clmark[secondLay][iSecondLay];
559                 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,mrk2);
560                 x2 = arr2->GetX();
561                 y2 = arr2->GetY();
562                 z2 = arr2->GetZ();
563                 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
564                 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
565                 phi2 = TMath::ATan2((y2-y1),(x2-x1));
566               } else { // special case of 1-point tracks, only for cosmics (B=0)
567                 x2 = primaryVertex[0];
568                 y2 = primaryVertex[1];
569                 z2 = primaryVertex[2];
570                 cv = 0;
571                 tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
572                 phi2 = TMath::ATan2((y1-y2),(x1-x2));
573               }
574
575               // create track and attach it the RecPoints
576               AliITStrackSA trac(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
577               for(Int_t iLay=5; iLay>=0; iLay--){
578                 Int_t iInLay=indices[iLay];
579                 AliITSRecPoint* cl=(AliITSRecPoint*)listlayer[iLay][iInLay];
580                 if(cl!=0){
581                   trac.AddClusterV2(iLay,(clind[iLay][iInLay] & 0x0fffffff)>>0);
582                   trac.AddClusterMark(iLay,clmark[iLay][iInLay]);
583                 }
584               }
585
586               //fit with Kalman filter using AliITStrackerMI::RefitAt()
587               AliITStrackSA ot(trac);
588
589               ot.ResetCovariance(10.);
590               ot.ResetClusters();
591               
592               // Propagate inside the innermost layer with a cluster 
593               if(ot.Propagate(ot.GetX()-0.1*ot.GetX())) {
594
595                 if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&ot,&trac)){ //fit from layer 1 to layer 6
596                   AliITStrackMI otrack2(ot);
597                   otrack2.ResetCovariance(10.); 
598                   otrack2.ResetClusters();
599                   //fit from layer 6 to layer 1
600                   if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),&otrack2,&ot)) {
601                     new(arrMI[nFoundTracks]) AliITStrackMI(otrack2);
602                     new(arrSA[nFoundTracks]) AliITStrackSA(trac);
603                     ++nFoundTracks;
604                   }
605                               
606                 }       
607               }
608             }//end loop layer 6
609           }//end loop layer 5
610         }//end loop layer 4        
611       }//end loop layer 3
612     }//end loop layer 2 
613   }//end loop layer 1
614
615
616
617
618   if(fListOfTracks->GetEntries()==0) return 0;
619
620   Int_t lowchi2 = FindTrackLowChiSquare();
621   AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
622   AliITStrackSA* trsa = (AliITStrackSA*)fListOfSATracks->At(lowchi2);
623  
624   if(otrack==0) return 0;
625
626   CookLabel(otrack,0.); //MI change - to see fake ratio
627   Int_t label=FindLabel(otrack);
628   otrack->SetLabel(label);  
629   Double_t low=0.;
630   Double_t up=0.51;    
631   otrack->CookdEdx(low,up);
632
633   //remove clusters of found track
634   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
635     for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
636       Int_t index = trsa->GetClusterMark(nlay,cln);
637       RemoveClusterCoord(nlay,index);
638     }    
639   }
640
641   return otrack;
642
643 }
644
645 //_______________________________________________________
646 void AliITStrackerSA::StoreTrack(AliITStrackV2 *t,AliESDEvent *event, Bool_t pureSA) const 
647 {
648   //
649   // Add new track to the ESD
650   //
651   AliESDtrack outtrack;
652   outtrack.UpdateTrackParams(t,AliESDtrack::kITSin);
653   if(pureSA) outtrack.SetStatus(AliESDtrack::kITSpureSA);
654   for(Int_t i=0;i<12;i++) {
655     outtrack.SetITSModuleIndex(i,t->GetModuleIndex(i));
656   }
657   Double_t sdedx[4]={0.,0.,0.,0.};
658   for(Int_t i=0; i<4; i++) sdedx[i]=t->GetSampledEdx(i);
659   outtrack.SetITSdEdxSamples(sdedx);
660
661
662   if(AliITSReconstructor::GetRecoParam()->GetSAUsedEdxInfo()){
663     Double_t mom=t->P();
664     Double_t ppid[AliPID::kSPECIES];
665     for(Int_t isp=0;isp<AliPID::kSPECIES;isp++) ppid[isp]=0.;
666     ppid[AliPID::kPion]=1.;
667     if(mom<0.7){
668       Double_t truncmean=t->GetdEdx();
669       Int_t ide=fITSPid->GetParticleIdFromdEdxVsP(mom,truncmean,kTRUE);
670       if(ide==AliPID::kProton){
671         ppid[AliPID::kProton]=1.;
672         ppid[AliPID::kPion]=0.;
673       }
674       else if(ide==AliPID::kKaon){ 
675         ppid[AliPID::kKaon]=1.; 
676         ppid[AliPID::kPion]=0.;
677       }
678     }
679     outtrack.SetITSpid(ppid);
680     outtrack.SetESDpid(ppid);    
681   }
682   event->AddTrack(&outtrack);
683
684   return;
685 }
686
687
688 //_______________________________________________________
689 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
690   //function used to to find the clusters associated to the track
691
692   if(ForceSkippingOfLayer(layer)) return 0;
693
694
695   Int_t nc=0;
696   AliITSlayer &lay = fgLayers[layer];
697   Double_t r=lay.GetR();
698   if(pflag==1){      
699     Double_t cx1,cx2,cy1,cy2;
700     FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
701     if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
702        return 0;
703     Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
704     Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
705     fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
706   }
707
708  
709   Double_t phiExpect=fPhiEstimate;
710   Double_t lamExpect=fLambdac;
711
712   Int_t ncl = fCluCoord[layer]->GetEntriesFast();
713   Int_t startcl=FindIndex(layer,lamExpect-lambdawindow*1.02);
714   Int_t endcl=FindIndex(layer,lamExpect+lambdawindow*1.02)+1;
715   if(endcl>=ncl) endcl=ncl-1;
716
717   for (Int_t index=startcl; index<=endcl; index++) {
718     //for (Int_t index=0; index<ncl; index++) {
719     AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
720
721
722     Double_t phi = arr->GetPhi();
723     Double_t deltaPhi = phi-phiExpect;
724     if(deltaPhi>TMath::Pi()) deltaPhi-=2*TMath::Pi();
725     else if(deltaPhi<-TMath::Pi()) deltaPhi+=2*TMath::Pi();
726     if (TMath::Abs(deltaPhi)>phiwindow) continue;
727     
728     Double_t lambda = arr->GetLambda();
729     if (TMath::Abs(lambda-lamExpect)>lambdawindow) continue;
730
731     if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
732     if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
733     Int_t orind = arr->GetOrInd();
734     trs->AddClusterSA(layer,orind);
735     trs->AddClusterMark(layer,index);
736     nc++;
737     fLambdac=lambda;
738     fPhiEstimate=phi;
739     
740     fPointc[0]=arr->GetX();
741     fPointc[1]=arr->GetY();
742     
743   }
744   return nc;
745 }
746
747 //________________________________________________________________
748 Bool_t AliITStrackerSA::SetFirstPoint(Int_t lay, Int_t clu, Double_t* primaryVertex){
749   // Sets the first point (seed) for tracking
750
751   AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay,clu);
752   fPhic = arr->GetPhi();
753   fLambdac = arr->GetLambda();
754   fPhiEstimate = fPhic;
755   fPoint1[0]=primaryVertex[0];
756   fPoint1[1]=primaryVertex[1];
757   fPoint2[0]=arr->GetX();
758   fPoint2[1]=arr->GetY();
759   return kTRUE; 
760 }
761
762 //________________________________________________________________
763 void AliITStrackerSA::UpdatePoints(){
764   //update of points for the estimation of the curvature  
765
766   fPoint2[0]=fPoint3[0];
767   fPoint2[1]=fPoint3[1];
768   fPoint3[0]=fPointc[0];
769   fPoint3[1]=fPointc[1];
770
771   
772 }
773
774 //___________________________________________________________________
775 Int_t AliITStrackerSA::FindEquation(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3,Double_t& a, Double_t& b, Double_t& c){
776
777    //given (x,y) of three recpoints (in global coordinates) 
778    //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
779
780    Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
781    if(den==0) return 0;
782    a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
783    b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
784    c = -x1*x1-y1*y1-a*x1-b*y1;
785    return 1;
786  }
787 //__________________________________________________________________________
788  Int_t AliITStrackerSA::FindIntersection(Double_t a1, Double_t b1, Double_t c1, Double_t c2,Double_t& x1,Double_t& y1, Double_t& x2, Double_t& y2){
789  
790  //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
791  //c2 is -rlayer*rlayer
792
793   if(a1==0) return 0;
794  Double_t m = c2-c1; 
795  Double_t aA = (b1*b1)/(a1*a1)+1;
796  Double_t bB = (-2*m*b1/(a1*a1));
797  Double_t cC = c2+(m*m)/(a1*a1);
798  Double_t dD = bB*bB-4*aA*cC;
799  if(dD<0) return 0;
800  
801  y1 = (-bB+TMath::Sqrt(dD))/(2*aA); 
802  y2 = (-bB-TMath::Sqrt(dD))/(2*aA); 
803  x1 = (c2-c1-b1*y1)/a1;
804  x2 = (c2-c1-b1*y2)/a1;
805
806  return 1; 
807 }
808 //____________________________________________________________________
809 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t 
810 x2,Double_t y2,Double_t x3,Double_t y3){
811
812   //calculates the curvature of track  
813   Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
814   if(den==0) return 0;
815   Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
816   Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
817   Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
818   Double_t xc=-a/2.;
819
820   if((a*a+b*b-4*c)<0) return 0;
821   Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
822   if(rad==0) return 0;
823   
824   if((x1>0 && y1>0 && x1<xc)) rad*=-1;
825   if((x1<0 && y1>0 && x1<xc)) rad*=-1;
826   //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
827   // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
828   
829   return 1/rad;
830  
831 }
832
833
834 //____________________________________________________________________
835 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
836
837   //Returns the point closest to pp
838
839   Double_t diff1 = p1-pp;
840   Double_t diff2 = p2-pp;
841   
842   if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
843   else fPhiEstimate=p2;  
844   return fPhiEstimate;
845   
846 }
847
848
849 //_________________________________________________________________
850 Int_t AliITStrackerSA::FindTrackLowChiSquare() const {
851   // returns track with lowest chi square  
852   Int_t dim=fListOfTracks->GetEntries();
853   if(dim<=1) return 0;
854   AliITStrackV2* trk = (AliITStrackV2*)fListOfTracks->At(0);
855   Double_t minChi2=trk->GetChi2();
856   Int_t index=0;
857   for(Int_t i=1;i<dim;i++){
858     trk = (AliITStrackV2*)fListOfTracks->At(i);
859     Double_t chi2=trk->GetChi2();
860     if(chi2<minChi2){
861       minChi2=chi2;
862       index=i;
863     }
864   }
865   return index;
866 }
867
868 //__________________________________________________________
869 Int_t AliITStrackerSA::FindLabel(AliITStrackV2* track){
870   // compute the track label starting from cluster labels
871   
872   Int_t labl[AliITSgeomTGeo::kNLayers][3];
873   Int_t cnts[AliITSgeomTGeo::kNLayers][3];
874   for(Int_t j=0;j<AliITSgeomTGeo::GetNLayers();j++){
875     for(Int_t k=0;k<3;k++){
876       labl[j][k]=-2;
877       cnts[j][k]=1;
878     }
879   }
880   Int_t iNotLabel=0;
881   for(Int_t i=0;i<track->GetNumberOfClusters(); i++) {
882     Int_t indexc = track->GetClusterIndex(i);
883     AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc);
884     Int_t iLayer=cl->GetLayer();
885     for(Int_t k=0;k<3;k++){
886       labl[iLayer][k]=cl->GetLabel(k);
887       if(labl[iLayer][k]<0) iNotLabel++;
888     }
889   }
890   if(iNotLabel==3*track->GetNumberOfClusters()) return -2;
891
892   for(Int_t j1=0;j1<AliITSgeomTGeo::kNLayers; j1++) {
893     for(Int_t j2=0; j2<j1;  j2++){
894       for(Int_t k1=0; k1<3; k1++){
895         for(Int_t k2=0; k2<3; k2++){
896           if(labl[j1][k1]>=0 && labl[j1][k1]==labl[j2][k2] && cnts[j2][k2]>0){
897             cnts[j2][k2]++;
898             cnts[j1][k1]=0;
899           }
900         }
901       }
902     }
903   }
904
905
906   Int_t cntMax=0;
907   Int_t label=-1;
908   for(Int_t j=0;j<AliITSgeomTGeo::kNLayers;j++){
909     for(Int_t k=0;k<3;k++){
910       if(cnts[j][k]>cntMax && labl[j][k]>=0){
911         cntMax=cnts[j][k];
912         label=labl[j][k];
913       }
914     }
915   }
916
917   Int_t lflag=0;
918   for(Int_t i=0;i<AliITSgeomTGeo::kNLayers;i++)
919     if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
920   
921   if(lflag<track->GetNumberOfClusters()) label = -label;
922   return label;
923 }
924 //_____________________________________________________________________________
925 void AliITStrackerSA::SetCalculatedWindowSizes(Int_t n, Double_t phimin, Double_t phimax, Double_t lambdamin, Double_t lambdamax){
926   // Set sizes of the phi and lambda windows used for track finding
927   fNloop = n;
928   if(fPhiWin) delete [] fPhiWin;
929   if(fLambdaWin) delete [] fLambdaWin;
930   fPhiWin = new Double_t[fNloop];
931   fLambdaWin = new Double_t[fNloop];
932   Double_t stepPhi=(phimax-phimin)/(Double_t)(fNloop-1);
933   Double_t stepLambda=(lambdamax-lambdamin)/(Double_t)(fNloop-1);
934   for(Int_t k=0;k<fNloop;k++){
935     Double_t phi=phimin+k*stepPhi;
936     Double_t lam=lambdamin+k*stepLambda;
937     fPhiWin[k]=phi;
938     fLambdaWin[k]=lam;
939   }
940 }
941 //_____________________________________________________________________________
942 void AliITStrackerSA::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
943   // Set sizes of the phi and lambda windows used for track finding
944   fNloop = n;
945   if(phi){ // user defined values
946     fPhiWin = new Double_t[fNloop];
947     fLambdaWin = new Double_t[fNloop];
948     for(Int_t k=0;k<fNloop;k++){
949       fPhiWin[k]=phi[k];
950       fLambdaWin[k]=lam[k];
951     }
952   }
953   else {  // default values
954             
955     Double_t phid[32]   = {0.002,0.003,0.004,0.0045,0.0047,
956                            0.005,0.0053,0.0055,0.006,0.0063,
957                            0.0065,0.007,0.0073,0.0075,0.0077,
958                            0.008,0.0083,0.0085,0.0087,0.009,
959                            0.0095,0.0097,0.01,0.0105,0.011,
960                            0.0115,0.012,0.0125,0.013,0.0135,
961                            0.0140,0.0145};
962     Double_t lambdad[32] = {0.003,0.004,0.005,0.005,0.005,
963                             0.005,0.005,0.006,0.006,0.006,
964                             0.006,0.007,0.007,0.007,0.007,
965                             0.007,0.007,0.007,0.007,0.007,
966                             0.007,0.007,0.008,0.008,0.008,
967                             0.008,0.008,0.008,0.008,0.008,
968                             0.008,0.008};
969     
970     if(fNloop!=32){
971       fNloop = 32;
972     }
973     
974     
975     fPhiWin = new Double_t[fNloop];
976     fLambdaWin = new Double_t[fNloop];
977
978     Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
979   
980     for(Int_t k=0;k<fNloop;k++){
981       fPhiWin[k]=phid[k]*factor;
982       fLambdaWin[k]=lambdad[k]*factor;
983     }
984   
985   }
986
987 }
988 //_______________________________________________________________________
989 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Double_t &x, Double_t &y,Double_t &z, const Double_t* vertex){
990   //Returns values of phi (azimuthal) and lambda angles for a given cluster
991 /*  
992   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
993   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
994   Double_t tx,ty,tz;  fGeom->GetTrans(lay,lad,det,tx,ty,tz);     
995
996   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
997   Double_t phi1=TMath::Pi()/2+alpha;
998   if (lay==1) phi1+=TMath::Pi();
999
1000   Double_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1001   Double_t r=tx*cp+ty*sp;
1002
1003   xyz= r*cp - cl->GetY()*sp;
1004   y= r*sp + cl->GetY()*cp;
1005   z=cl->GetZ();
1006 */
1007   Float_t xyz[3];
1008   cl->GetGlobalXYZ(xyz);
1009   x=xyz[0];
1010   y=xyz[1];
1011   z=xyz[2];
1012  
1013   phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1014   lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1015 }
1016
1017 //________________________________________________________________________
1018 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Double_t &sx,Double_t &sy, Double_t &sz){
1019
1020   //returns sigmax, y, z of cluster in global coordinates
1021 /*
1022   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1023   Int_t lay,lad,det; 
1024   AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1025  
1026   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1027   Double_t phi=TMath::Pi()/2+alpha;
1028   if (lay==1) phi+=TMath::Pi();
1029
1030   Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1031 */
1032   Float_t covm[6];
1033   cl->GetGlobalCov(covm);
1034   sx=TMath::Sqrt(covm[0]);
1035   sy=TMath::Sqrt(covm[3]);
1036   sz=TMath::Sqrt(covm[5]);
1037 /*
1038   sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1039   sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1040   sz = TMath::Sqrt(cl->GetSigmaZ2());
1041 */
1042 }
1043
1044 //________________________________________________________________________
1045 Int_t AliITStrackerSA::FindIndex(Int_t lay, Double_t lamVal) const {
1046   // Find the cluster at limit of lambda window 
1047
1048   Int_t base = 0;
1049   Int_t last = fCluCoord[lay]->GetEntriesFast()-1;
1050   if(last<0) return 0;
1051   Int_t position;
1052   Double_t lamfirst=((AliITSclusterTable*)fCluCoord[lay]->At(base))->GetLambda();
1053   if(lamfirst>lamVal) return base;
1054   Double_t lamlast=((AliITSclusterTable*)fCluCoord[lay]->At(last))->GetLambda();
1055   if(lamlast<=lamVal) return last;
1056   while (last >= base) {
1057     position = (base+last) / 2;
1058     Double_t a=((AliITSclusterTable*)fCluCoord[lay]->At(position))->GetLambda()-lamVal;
1059     Double_t b=((AliITSclusterTable*)fCluCoord[lay]->At(position+1))->GetLambda()-lamVal;
1060     if(a*b<=0) return position;
1061     if(a>0) last = position;
1062     else  base = position;
1063   }
1064   return 0;
1065 }
1066