Updated code for tracking V2 (from Y. Belikov)
[u/mrichter/AliRoot.git] / ITS / AliITStrackerV2.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 //-------------------------------------------------------------------------
17 //               Implementation of the ITS tracker class
18 //
19 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
20 //     dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
21 //-------------------------------------------------------------------------
22 #include <TFile.h>
23 #include <TTree.h>
24 #include <TRandom.h>
25 #include <iostream.h>
26
27 #include "AliITSgeom.h"
28 #include "AliITSRecPoint.h"
29 #include "AliTPCtrack.h"
30 #include "AliITSclusterV2.h"
31 #include "AliITStrackerV2.h"
32
33 //#define DEBUG
34
35 #ifdef DEBUG
36 Int_t LAB=5126;
37 #endif
38
39 AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers
40
41 AliITStrackerV2::
42 AliITStrackerV2(const AliITSgeom *geom, Int_t eventn) throw (const Char_t *) :
43 AliTracker() {
44   //--------------------------------------------------------------------
45   //This is the AliITStracker constructor
46   //It also reads clusters from a root file
47   //--------------------------------------------------------------------
48   fEventN=eventn;  //MI change add event number  - used to generate identifier 
49
50   AliITSgeom *g=(AliITSgeom*)geom;
51
52   Float_t x,y,z;
53   Int_t i;
54   for (i=1; i<kMaxLayer+1; i++) {
55     Int_t nlad=g->GetNladders(i);
56     Int_t ndet=g->GetNdetectors(i);
57
58     g->GetTrans(i,1,1,x,y,z); 
59     Double_t r=TMath::Sqrt(x*x + y*y);
60     Double_t poff=TMath::ATan2(y,x);
61     Double_t zoff=z;
62
63     g->GetTrans(i,1,2,x,y,z);
64     r += TMath::Sqrt(x*x + y*y);
65     g->GetTrans(i,2,1,x,y,z);
66     r += TMath::Sqrt(x*x + y*y);
67     g->GetTrans(i,2,2,x,y,z);
68     r += TMath::Sqrt(x*x + y*y);
69     r*=0.25;
70
71     new (fLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
72
73     for (Int_t j=1; j<nlad+1; j++) {
74       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
75         Float_t x,y,zshift; g->GetTrans(i,j,k,x,y,zshift); 
76         Double_t rot[9]; g->GetRotMatrix(i,j,k,rot);
77
78         Double_t r =-x*rot[1] + y*rot[0];         if (i==1) r=-r;
79         Double_t phi=TMath::ATan2(rot[1],rot[0]); if (i==1) phi-=3.1415927;
80         phi+=0.5*TMath::Pi(); 
81         AliITSdetector &det=fLayers[i-1].GetDetector((j-1)*ndet + k-1); 
82
83 if (phi<0) phi += 2*TMath::Pi();
84
85         new(&det) AliITSdetector(r,phi); 
86       } 
87     }  
88
89   }
90
91   try {
92      //Read clusters
93      //MI change 
94      char   cname[100]; 
95      sprintf(cname,"TreeC_ITS_%d",eventn);
96      TTree *cTree=(TTree*)gDirectory->Get(cname);
97
98      if (!cTree) throw 
99         ("AliITStrackerV2::AliITStrackerV2 can't get cTree !\n");
100
101      TBranch *branch=cTree->GetBranch("Clusters");
102      if (!branch) throw 
103         ("AliITStrackerV2::AliITStrackerV2 can't get Clusters branch !\n");
104
105      TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
106      branch->SetAddress(&clusters);
107
108      Int_t nentr=(Int_t)cTree->GetEntries();
109      for (i=0; i<nentr; i++) {
110        if (!cTree->GetEvent(i)) continue;
111        Int_t lay,lad,det; g->GetModuleId(i,lay,lad,det);
112        Int_t ncl=clusters->GetEntriesFast();
113        while (ncl--) {
114          AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
115
116 #ifdef DEBUG
117 if (c->GetLabel(2)!=LAB)
118 if (c->GetLabel(1)!=LAB)
119 if (c->GetLabel(0)!=LAB) continue;
120 cout<<lay-1<<' '<<lad-1<<' '<<det-1<<' '<<c->GetY()<<' '<<c->GetZ()<<endl;
121 #endif
122
123          fLayers[lay-1].InsertCluster(new AliITSclusterV2(*c));
124        }
125        clusters->Delete();
126      }
127      delete cTree; //Thanks to Mariana Bondila
128   }
129   catch (const Char_t *msg) {
130     cerr<<msg<<endl;
131     throw;
132   }
133
134   fI=kMaxLayer;
135
136   fPass=0;
137   fConstraint[0]=1; fConstraint[1]=0;
138 }
139
140 #ifdef DEBUG
141 static Int_t lbl;
142 #endif
143
144 Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) {
145   //--------------------------------------------------------------------
146   //This functions reconstructs ITS tracks
147   //--------------------------------------------------------------------
148   TFile *in=(TFile*)inp;
149   TDirectory *savedir=gDirectory; 
150
151   if (!in->IsOpen()) {
152      cerr<<"AliITStrackerV2::Clusters2Tracks(): ";
153      cerr<<"file with TPC tracks is not open !\n";
154      return 1;
155   }
156
157   if (!out->IsOpen()) {
158      cerr<<"AliITStrackerV2::Clusters2Tracks(): ";
159      cerr<<"file for ITS tracks is not open !\n";
160      return 2;
161   }
162
163   in->cd();
164  
165   char   tname[100];
166   Int_t nentr=0; TObjArray itsTracks(10000);
167
168   {/* Read TPC tracks */ 
169     sprintf(tname,"TreeT_TPC_%d",fEventN);
170     TTree *tpcTree=(TTree*)in->Get(tname);
171     if (!tpcTree) {
172        cerr<<"AliITStrackerV2::Clusters2Tracks(): "
173              "can't get a tree with TPC tracks !\n";
174        return 3;
175     }
176     AliTPCtrack *itrack=new AliTPCtrack; 
177     tpcTree->SetBranchAddress("tracks",&itrack);
178     nentr=(Int_t)tpcTree->GetEntries();
179     for (Int_t i=0; i<nentr; i++) {
180        tpcTree->GetEvent(i);
181        AliITStrackV2 *t=new AliITStrackV2(*itrack);
182        if (TMath::Abs(t->GetD())>4) continue;
183
184        t->PropagateTo(80.,0.0053);
185        if (TMath::Abs(t->GetY())>13.) t->CorrectForMaterial(0.03);
186        if (TMath::Abs(t->GetZ())<0.2) t->CorrectForMaterial(0.40);
187        t->PropagateTo(61.,0.0052);
188        Double_t xk=52.,x,y,z; t->GetGlobalXYZat(xk,x,y,z);
189        if (TMath::Abs(y)<7.77) t->PropagateTo(xk,0.19,24.); 
190        t->PropagateTo(50.,0.001);
191
192        itsTracks.AddLast(t);
193     }
194     delete tpcTree; //Thanks to Mariana Bondila
195     delete itrack;
196   }
197   itsTracks.Sort();
198
199   out->cd();
200
201   sprintf(tname,"TreeT_ITS_%d",fEventN);
202   TTree itsTree(tname,"Tree with ITS tracks");
203   AliITStrackV2 *otrack=&fBestTrack;
204   itsTree.Branch("tracks","AliITStrackV2",&otrack,32000,0);
205
206   for (fPass=0; fPass<2; fPass++) {
207      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
208      for (Int_t i=0; i<nentr; i++) {
209        if (i%10==0) cerr<<nentr-i<<" \r";
210        AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i);
211        if (t==0) continue;           //this track has been already tracked
212        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
213
214 #ifdef DEBUG
215 lbl=tpcLabel;
216 if (TMath::Abs(tpcLabel)!=LAB) continue;
217 cout<<tpcLabel<<" *****************\n";
218 #endif
219
220        try {
221          ResetTrackToFollow(*t);
222        } catch (const Char_t *msg) {
223          cerr<<msg<<endl;
224          continue;
225        }
226        ResetBestTrack();
227
228        for (FollowProlongation(); fI<kMaxLayer; fI++) {
229           while (TakeNextProlongation()) FollowProlongation();
230        }
231
232 #ifdef DEBUG
233 cout<<fBestTrack.GetNumberOfClusters()<<" number of clusters\n\n";
234 #endif
235
236        if (fBestTrack.GetNumberOfClusters() >= kMaxLayer-kLayersToSkip) {
237           fBestTrack.SetLabel(tpcLabel);
238           fBestTrack.CookdEdx();
239           CookLabel(&fBestTrack,0.); //For comparison only
240           itsTree.Fill();
241           UseClusters(&fBestTrack);
242           delete itsTracks.RemoveAt(i);
243        }
244
245      }
246   }
247
248   itsTree.Write();
249
250   itsTracks.Delete();
251   savedir->cd();
252   cerr<<"Number of TPC tracks: "<<nentr<<endl;
253   cerr<<"Number of prolonged tracks: "<<itsTree.GetEntries()<<endl;
254
255   return 0;
256 }
257
258 Int_t AliITStrackerV2::PropagateBack(const TFile *inp, TFile *out) {
259   //--------------------------------------------------------------------
260   //This functions propagates reconstructed ITS tracks back
261   //--------------------------------------------------------------------
262   TFile *in=(TFile*)inp;
263   TDirectory *savedir=gDirectory; 
264
265   if (!in->IsOpen()) {
266      cerr<<"AliITStrackerV2::PropagateBack(): ";
267      cerr<<"file with ITS tracks is not open !\n";
268      return 1;
269   }
270
271   if (!out->IsOpen()) {
272      cerr<<"AliITStrackerV2::PropagateBack(): ";
273      cerr<<"file for back propagated ITS tracks is not open !\n";
274      return 2;
275   }
276
277   in->cd();
278   TTree *itsTree=(TTree*)in->Get("TreeT_ITS_0");
279   if (!itsTree) {
280      cerr<<"AliITStrackerV2::PropagateBack() ";
281      cerr<<"can't get a tree with ITS tracks !\n";
282      return 3;
283   }
284   AliITStrackV2 *itrack=new AliITStrackV2; 
285   itsTree->SetBranchAddress("tracks",&itrack);
286
287   out->cd();
288   TTree backTree("TreeT_ITSb_0","Tree with back propagated ITS tracks");
289   AliTPCtrack *otrack=0;
290   backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
291
292   Int_t ntrk=0;
293
294   Int_t nentr=(Int_t)itsTree->GetEntries();
295   for (Int_t i=0; i<nentr; i++) {
296     itsTree->GetEvent(i);
297     ResetTrackToFollow(*itrack);
298     fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters();
299     Int_t itsLabel=fTrackToFollow.GetLabel(); //save the ITS track label
300
301 #ifdef DEBUG
302 if (TMath::Abs(itsLabel)!=LAB) continue;
303 cout<<itsLabel<<" *****************\n";
304 #endif
305
306     try {
307        Int_t nc=itrack->GetNumberOfClusters();
308 #ifdef DEBUG
309 for (Int_t k=0; k<nc; k++) {
310     Int_t index=itrack->GetClusterIndex(k);
311     AliITSclusterV2 *c=(AliITSclusterV2*)GetCluster(index);
312     cout<<c->GetLabel(0)<<' '<<c->GetLabel(1)<<' '<<c->GetLabel(2)<<endl;
313 }
314 #endif       
315        Int_t idx=0, l=0; 
316        const  AliITSclusterV2 *c=0; 
317        if (nc--) {
318           idx=itrack->GetClusterIndex(nc); l=(idx&0xf0000000)>>28;
319           c=(AliITSclusterV2*)GetCluster(idx);
320        }
321        for (fI=0; fI<kMaxLayer; fI++) {
322          AliITSlayer &layer=fLayers[fI];
323          Double_t r=layer.GetR();
324          if (fI==2 || fI==4) {             
325             Double_t rs=0.5*(fLayers[fI-1].GetR() + r);
326             Double_t ds=0.034; if (fI==4) ds=0.039;
327             Double_t dx0r=ds/21.82*2.33, dr=ds*2.33;
328             if (!fTrackToFollow.PropagateTo(rs,1*dx0r,dr)) throw "";
329          }
330
331          Double_t x,y,z;
332          if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) 
333             throw "AliITStrackerV2::PropagateBack: failed to estimate track !";
334          Double_t phi=TMath::ATan2(y,x);
335          Int_t idet=layer.FindDetectorIndex(phi,z);
336          if (idet<0) 
337          throw "AliITStrackerV2::PropagateBack: failed to find a detector !\n";
338          const AliITSdetector &det=layer.GetDetector(idet);
339          r=det.GetR(); phi=det.GetPhi();
340          if (!fTrackToFollow.Propagate(phi,r)) throw "";
341          fTrackToFollow.SetDetectorIndex(idet);
342
343          const AliITSclusterV2 *cl=0;
344          Int_t index=0;
345          Double_t maxchi2=kMaxChi2;
346
347          if (l==fI) {
348            idet=c->GetDetectorIndex();
349            if (idet != fTrackToFollow.GetDetectorIndex()) {
350              const AliITSdetector &det=layer.GetDetector(idet);
351              r=det.GetR(); phi=det.GetPhi();
352              if (!fTrackToFollow.Propagate(phi,r)) throw "";
353              fTrackToFollow.SetDetectorIndex(idet);
354            }
355            Double_t chi2=fTrackToFollow.GetPredictedChi2(c);
356            if (chi2<kMaxChi2) {
357               cl=c; maxchi2=chi2; index=idx;
358            }
359            if (nc--) {
360               idx=itrack->GetClusterIndex(nc); l=(idx&0xf0000000)>>28;
361               c=(AliITSclusterV2*)GetCluster(idx);
362            } 
363          }
364
365          if (fTrackToFollow.GetNumberOfClusters()>2) {
366            Double_t dz=3*TMath::Sqrt(fTrackToFollow.GetSigmaZ2()+kSigmaZ2[fI]);
367            Double_t dy=3*TMath::Sqrt(fTrackToFollow.GetSigmaY2()+kSigmaY2[fI]);
368            Double_t zmin=fTrackToFollow.GetZ() - dz;
369            Double_t zmax=fTrackToFollow.GetZ() + dz;
370            Double_t ymin=fTrackToFollow.GetY() + phi*r - dy;
371            Double_t ymax=fTrackToFollow.GetY() + phi*r + dy;
372            layer.SelectClusters(zmin,zmax,ymin,ymax);
373
374            const AliITSclusterV2 *cc=0; Int_t ci;
375            while ((cc=layer.GetNextCluster(ci))!=0) {
376               idet=cc->GetDetectorIndex();
377               if (idet != fTrackToFollow.GetDetectorIndex()) continue;
378               Double_t chi2=fTrackToFollow.GetPredictedChi2(cc);
379               if (chi2<maxchi2) {
380                  cl=cc; index=(fI<<28)+ci; maxchi2=chi2;
381               }
382            }
383          }
384
385          if (cl) {
386             if (!fTrackToFollow.Update(cl,maxchi2,index)) 
387               cerr<<"AliITStrackerV2::PropagateBack: filtering failed !\n";
388          }
389
390          x=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ());
391          fTrackToFollow.CorrectForMaterial(-x); 
392
393        }
394
395        Double_t xk=52.,x,y,z; fTrackToFollow.GetGlobalXYZat(xk,x,y,z);
396        if (TMath::Abs(y)<7.77) 
397           fTrackToFollow.PropagateTo(xk,-0.19,24.); 
398        fTrackToFollow.PropagateTo(61,-0.0110);
399        fTrackToFollow.PropagateTo(80.,-0.0053);
400
401        fTrackToFollow.SetLabel(itsLabel);
402        otrack=new AliTPCtrack(fTrackToFollow,fTrackToFollow.GetAlpha()); 
403        backTree.Fill(); delete otrack;
404        UseClusters(&fTrackToFollow);
405        cerr << ++ntrk << "                \r";
406     }
407     catch (const Char_t *msg) {
408        cerr<<msg<<endl;
409     }
410   }
411
412   backTree.Write();
413   savedir->cd();
414   cerr<<"Number of ITS tracks: "<<nentr<<endl;
415   cerr<<"Number of back propagated ITS tracks: "<<ntrk<<endl;
416
417   delete itrack;
418
419   delete itsTree; //Thanks to Mariana Bondila
420
421   return 0;
422 }
423
424
425 AliCluster *AliITStrackerV2::GetCluster(Int_t index) const {
426   //--------------------------------------------------------------------
427   //       Return pointer to a given cluster
428   //--------------------------------------------------------------------
429   Int_t l=(index & 0xf0000000) >> 28;
430   Int_t c=(index & 0x0fffffff) >> 00;
431   return fLayers[l].GetCluster(c);
432 }
433
434
435 void AliITStrackerV2::FollowProlongation() {
436   //--------------------------------------------------------------------
437   //This function finds a track prolongation 
438   //--------------------------------------------------------------------
439   Int_t tryAgain=kLayersToSkip;
440
441   while (fI) {
442     Int_t i=fI-1;
443 #ifdef DEBUG
444 cout<<i<<' ';
445 #endif
446     AliITSlayer &layer=fLayers[i];
447     AliITStrackV2 &track=fTracks[i];
448
449     Double_t r=layer.GetR();
450     if (i==3 || i==1) {
451        Double_t rs=0.5*(fLayers[i+1].GetR() + r);
452        Double_t d=0.011; if (i==3) d=0.0053;
453        if (!fTrackToFollow.PropagateTo(rs,d)) {
454          //cerr<<"AliITStrackerV2::FollowProlongation: "
455          //"propagation failed !\n";
456          break;
457        }
458     }
459
460     //find intersection
461     Double_t x,y,z;  
462     if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) {
463       //cerr<<"AliITStrackerV2::FollowProlongation: "
464       //"failed to estimate track !\n";
465       break;
466     }
467     Double_t phi=TMath::ATan2(y,x);
468     Int_t idet=layer.FindDetectorIndex(phi,z);
469     if (idet<0) {
470       //cerr<<"AliITStrackerV2::FollowProlongation: "
471       //"failed to find a detector !\n";
472       break;
473     }
474
475     //propagate to the intersection
476     const AliITSdetector &det=layer.GetDetector(idet);
477     phi=det.GetPhi();
478     if (!fTrackToFollow.Propagate(phi,det.GetR())) {
479       //cerr<<"AliITStrackerV2::FollowProlongation: propagation failed !\n";
480       break;
481     }
482     fTrackToFollow.SetDetectorIndex(idet);
483
484     //Select possible prolongations and store the current track estimation
485     track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow);
486     Double_t dz=3*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]);
487     if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl());
488     if (dz > kMaxRoad) {
489       //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Z !\n";
490       break;
491     }
492
493     if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+dz) break;
494
495     Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]);
496     if (dy < 0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp());
497     if (dy > kMaxRoad) {
498       //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Y !\n";
499       break;
500     }
501
502     Double_t zmin=track.GetZ() - dz; 
503     Double_t zmax=track.GetZ() + dz;
504     Double_t ymin=track.GetY() + r*phi - dy;
505     Double_t ymax=track.GetY() + r*phi + dy;
506     layer.SelectClusters(zmin,zmax,ymin,ymax); 
507     fI--;
508
509     //take another prolongation
510     if (!TakeNextProlongation()) if (!tryAgain--) break;
511     tryAgain=kLayersToSkip;
512
513   } 
514
515   //deal with the best track
516   Int_t ncl=fTrackToFollow.GetNumberOfClusters();
517   Int_t nclb=fBestTrack.GetNumberOfClusters();
518   if (ncl)
519   if (ncl >= nclb) {
520      Double_t chi2=fTrackToFollow.GetChi2();
521      if (chi2/ncl < kChi2PerCluster) {        
522         if (ncl > nclb || chi2 < fBestTrack.GetChi2()) {
523            ResetBestTrack();
524         }
525      }
526   }
527
528 }
529
530 Int_t AliITStrackerV2::TakeNextProlongation() {
531   //--------------------------------------------------------------------
532   // This function takes another track prolongation 
533   //
534   //  dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch 
535   //--------------------------------------------------------------------
536   AliITSlayer &layer=fLayers[fI];
537   AliITStrackV2 &t=fTracks[fI];
538
539   Double_t dz=4*TMath::Sqrt(t.GetSigmaZ2() + kSigmaZ2[fI]);
540   Double_t dy=4*TMath::Sqrt(t.GetSigmaY2() + kSigmaY2[fI]);
541
542   const AliITSclusterV2 *c=0; Int_t ci=-1;
543   Double_t chi2=12345.;
544   while ((c=layer.GetNextCluster(ci))!=0) {
545     //if (c->GetLabel(0)!=TMath::Abs(lbl)) continue; 
546     Int_t idet=c->GetDetectorIndex();
547
548     if (t.GetDetectorIndex()!=idet) {
549        const AliITSdetector &det=layer.GetDetector(idet);
550        if (!t.Propagate(det.GetPhi(),det.GetR())) {
551          //cerr<<"AliITStrackerV2::TakeNextProlongation: "
552          //"propagation failed !\n";
553          continue;
554        }
555        t.SetDetectorIndex(idet);
556        if (TMath::Abs(t.GetZ()-GetZ()) > layer.GetR()+dz) continue;
557
558 #ifdef DEBUG
559 cout<<fI<<" change detector !\n";
560 #endif
561
562     }
563
564     if (TMath::Abs(t.GetZ() - c->GetZ()) > dz) continue;
565     if (TMath::Abs(t.GetY() - c->GetY()) > dy) continue;
566
567     chi2=t.GetPredictedChi2(c); if (chi2<kMaxChi2) break;
568   }
569
570 #ifdef DEBUG
571 cout<<fI<<" chi2="<<chi2<<' '<<t.GetY()<<' '<<t.GetZ()<<' '<<dy<<' '<<dz<<endl;
572 #endif
573
574   if (chi2>=kMaxChi2) return 0;
575   if (!c) return 0;
576
577   ResetTrackToFollow(t);
578
579   if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) {
580      //cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
581      return 0;
582   }
583
584   fTrackToFollow.
585     SetSampledEdx(c->GetQ(),fTrackToFollow.GetNumberOfClusters()-1); //b.b.
586
587   {
588    Double_t d=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ());
589    fTrackToFollow.CorrectForMaterial(d);
590   }
591   if (fConstraint[fPass]) {
592      Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
593      fTrackToFollow.Improve(d,GetY(),GetZ());
594   }
595
596 #ifdef DEBUG
597 cout<<"accepted lab="<<c->GetLabel(0)<<' '
598     <<fTrackToFollow.GetNumberOfClusters()<<' '
599     <<fTrackToFollow.GetY()<<' '<<fTrackToFollow.GetZ()<<' '
600     <<fTrackToFollow.Get1Pt()<<endl<<endl;
601 #endif
602
603   return 1;
604 }
605
606
607 AliITStrackerV2::AliITSlayer::AliITSlayer() {
608   //--------------------------------------------------------------------
609   //default AliITSlayer constructor
610   //--------------------------------------------------------------------
611   fN=0;
612   fDetectors=0;
613 }
614
615 AliITStrackerV2::AliITSlayer::
616 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd) {
617   //--------------------------------------------------------------------
618   //main AliITSlayer constructor
619   //--------------------------------------------------------------------
620   fR=r; fPhiOffset=p; fZOffset=z;
621   fNladders=nl; fNdetectors=nd;
622   fDetectors=new AliITSdetector[fNladders*fNdetectors];
623
624   fN=0;
625   fI=0;
626 }
627
628 AliITStrackerV2::AliITSlayer::~AliITSlayer() {
629   //--------------------------------------------------------------------
630   // AliITSlayer destructor
631   //--------------------------------------------------------------------
632   delete[] fDetectors;
633   for (Int_t i=0; i<fN; i++) delete fClusters[i];
634 }
635
636 Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSclusterV2 *c) {
637   //--------------------------------------------------------------------
638   //This function adds a cluster to this layer
639   //--------------------------------------------------------------------
640   if (fN==kMaxClusterPerLayer) {
641      cerr<<"AliITStrackerV2::AliITSlayer::InsertCluster(): "
642            "Too many clusters !\n";
643      return 1;
644   }
645
646   if (fN==0) {fClusters[fN++]=c; return 0;}
647   Int_t i=FindClusterIndex(c->GetZ());
648   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliITSclusterV2*));
649   fClusters[i]=c; fN++;
650
651   return 0;
652 }
653
654 Int_t AliITStrackerV2::AliITSlayer::FindClusterIndex(Double_t z) const {
655   //--------------------------------------------------------------------
656   // This function returns the index of the nearest cluster 
657   //--------------------------------------------------------------------
658   if (fN==0) return 0;
659   if (z <= fClusters[0]->GetZ()) return 0;
660   if (z > fClusters[fN-1]->GetZ()) return fN;
661   Int_t b=0, e=fN-1, m=(b+e)/2;
662   for (; b<e; m=(b+e)/2) {
663     if (z > fClusters[m]->GetZ()) b=m+1;
664     else e=m; 
665   }
666   return m;
667 }
668
669 void AliITStrackerV2::AliITSlayer::
670 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
671   //--------------------------------------------------------------------
672   // This function sets the "window"
673   //--------------------------------------------------------------------
674   fI=FindClusterIndex(zmin); fZmax=zmax;
675   Double_t circle=2*TMath::Pi()*fR;
676   if (ymax>circle) { ymax-=circle; ymin-=circle; }
677   fYmin=ymin; fYmax=ymax;
678 }
679
680 const AliITSclusterV2 *AliITStrackerV2::AliITSlayer::GetNextCluster(Int_t &ci){
681   //--------------------------------------------------------------------
682   // This function returns clusters within the "window" 
683   //--------------------------------------------------------------------
684   const AliITSclusterV2 *cluster=0;
685   for (Int_t i=fI; i<fN; i++) {
686     const AliITSclusterV2 *c=fClusters[i];
687     if (c->GetZ() > fZmax) break;
688     if (c->IsUsed()) continue;
689     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());    
690     Double_t y=fR*det.GetPhi() + c->GetY();
691
692     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
693     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
694
695     if (y<fYmin) continue;
696     if (y>fYmax) continue;
697     cluster=c; ci=i;
698     fI=i+1;
699     break; 
700   }
701
702   return cluster;
703 }
704
705 Int_t AliITStrackerV2::AliITSlayer::
706 FindDetectorIndex(Double_t phi, Double_t z) const {
707   //--------------------------------------------------------------------
708   //This function finds the detector crossed by the track
709   //--------------------------------------------------------------------
710   Double_t dphi=phi-fPhiOffset;
711   if      (dphi <  0) dphi += 2*TMath::Pi();
712   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
713   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
714
715   Double_t dz=fZOffset-z;
716   Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
717
718   if (np>=fNladders) np-=fNladders;
719   if (np<0)          np+=fNladders;
720
721 #ifdef DEBUG
722 cout<<np<<' '<<nz<<endl;
723 #endif
724
725   return np*fNdetectors + nz;
726 }
727
728 Double_t 
729 AliITStrackerV2::AliITSlayer::GetThickness(Double_t y, Double_t z) const {
730   //--------------------------------------------------------------------
731   //This function returns the layer thickness at this point (units X0)
732   //--------------------------------------------------------------------
733   Double_t d=0.0085;
734
735   if (43<fR&&fR<45) { //SSD2
736      d=0.0036;
737      if (TMath::Abs(y-0.00)>3.40) d+=0.0036;
738      if (TMath::Abs(y-2.50)<0.10) d+=(0.02-0.0036);
739      if (TMath::Abs(y+1.90)<0.10) d+=(0.02-0.0036);
740      for (Int_t i=0; i<12; i++) {
741        if (TMath::Abs(z-3.6*(i+0.5))<0.20) {d+=0.0036; break;}
742        if (TMath::Abs(z+3.6*(i+0.5))<0.20) {d+=0.0036; break;}         
743        if (TMath::Abs(z-3.6*(i+0.929))<0.50) {d+=(0.02-0.0036); break;}
744        if (TMath::Abs(z+3.6*(i+0.104))<0.50) {d+=(0.02-0.0036); break;}
745      }
746   } else 
747   if (37<fR&&fR<41) { //SSD1
748      d=0.0036;
749      if (TMath::Abs(y-0.00)>3.40) d+=0.0036;
750      if (TMath::Abs(y-2.20)<0.10) d+=(0.02-0.0036);
751      if (TMath::Abs(y+2.20)<0.10) d+=(0.02-0.0036);
752      for (Int_t i=0; i<11; i++) {
753        if (TMath::Abs(z-3.6*i)<0.20) {d+=0.0036; break;}
754        if (TMath::Abs(z+3.6*i)<0.20) {d+=0.0036; break;}         
755        if (TMath::Abs(z-3.6*(i+0.54))<0.50) {d+=(0.02-0.0036); break;}
756        if (TMath::Abs(z+3.6*(i+0.58))<0.50) {d+=(0.02-0.0036); break;}         
757      }
758   } else
759   if (13<fR&&fR<26) { //SDD
760      d=0.0034;
761      if (TMath::Abs(y-0.00)>3.30) d+=0.0034;
762      if (TMath::Abs(y-2.10)<0.20) d+=0.0034*3;
763      if (TMath::Abs(y+2.10)<0.20) d+=0.0034*3;
764      for (Int_t i=0; i<4; i++) { 
765        if (TMath::Abs(z-7.3*i)<0.60) {d+=0.0034; break;}
766        if (TMath::Abs(z+7.3*i)<0.60) {d+=0.0034; break;}
767      }
768   } else
769   if (6<fR&&fR<8) {   //SPD2
770      d=0.0093;
771      if (TMath::Abs(y-3.08)>0.45) d+=0.0064;
772      if (TMath::Abs(y-3.03)<0.10) d+=0.0192;
773   } else
774   if (3<fR&&fR<5) {   //SPD1
775      d=0.0093;
776      if (TMath::Abs(y+0.21)>0.55) d+=0.0064;
777      if (TMath::Abs(y+0.10)<0.10) d+=0.0192;
778   }
779
780   d+=0.002;
781
782   return d;
783 }
784
785 Double_t AliITStrackerV2::GetEffectiveThickness(Double_t y,Double_t z) const
786 {
787   //--------------------------------------------------------------------
788   //Returns the thickness between the current layer and the vertex (units X0)
789   //--------------------------------------------------------------------
790   Double_t d=0.1/65.19*1.848;
791
792   Double_t xn=fLayers[fI].GetR();
793   for (Int_t i=0; i<fI; i++) {
794     Double_t xi=fLayers[i].GetR();
795     d+=fLayers[i].GetThickness(y,z)*xi*xi;
796   }
797
798   if (fI>1) {
799     Double_t xi=0.5*(fLayers[1].GetR()+fLayers[2].GetR());
800     d+=0.011*xi*xi;
801   }
802
803   if (fI>3) {
804     Double_t xi=0.5*(fLayers[3].GetR()+fLayers[4].GetR());
805     d+=0.0053*xi*xi;
806   }
807   return d/(xn*xn);
808 }
809
810
811
812 Int_t AliITStrackerV2::AliITSlayer::InRoad() const {
813   //--------------------------------------------------------------------
814   // This function returns number of clusters within the "window" 
815   //--------------------------------------------------------------------
816   Int_t ncl=0;
817   for (Int_t i=fI; i<fN; i++) {
818     const AliITSclusterV2 *c=fClusters[i];
819     if (c->GetZ() > fZmax) break;
820     if (c->IsUsed()) continue;
821     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());    
822     Double_t y=fR*det.GetPhi() + c->GetY();
823
824     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
825     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
826
827     if (y<fYmin) continue;
828     if (y>fYmax) continue;
829     ncl++;
830   }
831   return ncl;
832 }
833
834
835