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