8c6314b9d7a21aa3c181b036289300cc99732f85
[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 <Riostream.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=113;
37 #endif
38
39 AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers
40
41 AliITStrackerV2::AliITStrackerV2(const AliITSgeom *geom) : AliTracker() {
42   //--------------------------------------------------------------------
43   //This is the AliITStrackerV2 constructor
44   //--------------------------------------------------------------------
45   AliITSgeom *g=(AliITSgeom*)geom;
46
47   Float_t x,y,z;
48   Int_t i;
49   for (i=1; i<kMaxLayer+1; i++) {
50     Int_t nlad=g->GetNladders(i);
51     Int_t ndet=g->GetNdetectors(i);
52
53     g->GetTrans(i,1,1,x,y,z); 
54     Double_t r=TMath::Sqrt(x*x + y*y);
55     Double_t poff=TMath::ATan2(y,x);
56     Double_t zoff=z;
57
58     g->GetTrans(i,1,2,x,y,z);
59     r += TMath::Sqrt(x*x + y*y);
60     g->GetTrans(i,2,1,x,y,z);
61     r += TMath::Sqrt(x*x + y*y);
62     g->GetTrans(i,2,2,x,y,z);
63     r += TMath::Sqrt(x*x + y*y);
64     r*=0.25;
65
66     new (fLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
67
68     for (Int_t j=1; j<nlad+1; j++) {
69       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
70         Float_t x,y,zshift; g->GetTrans(i,j,k,x,y,zshift); 
71         Double_t rot[9]; g->GetRotMatrix(i,j,k,rot);
72
73         Double_t r =-x*rot[1] + y*rot[0];         if (i==1) r=-r;
74         Double_t phi=TMath::ATan2(rot[1],rot[0]); if (i==1) phi-=3.1415927;
75         phi+=0.5*TMath::Pi(); if (phi<0) phi += 2*TMath::Pi();
76         AliITSdetector &det=fLayers[i-1].GetDetector((j-1)*ndet + k-1); 
77
78         new(&det) AliITSdetector(r,phi); 
79       } 
80     }  
81
82   }
83
84   fI=kMaxLayer;
85
86   fPass=0;
87   fConstraint[0]=1; fConstraint[1]=0;
88 }
89
90 void AliITStrackerV2::LoadClusters() {
91   //--------------------------------------------------------------------
92   //This function loads ITS clusters
93   //--------------------------------------------------------------------
94   char   cname[100]; 
95   sprintf(cname,"TreeC_ITS_%d",GetEventNumber());
96   TTree *cTree=(TTree*)gDirectory->Get(cname);
97   if (!cTree) { 
98     cerr<<"AliITStrackerV2::LoadClusters can't get cTree !\n";
99     exit(1);
100   }
101   TBranch *branch=cTree->GetBranch("Clusters");
102   if (!branch) { 
103     cerr<<"AliITStrackerV2::LoadClusters can't get Clusters branch !\n";
104     exit(1);
105   }
106
107   TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
108   branch->SetAddress(&clusters);
109
110   Int_t j=0;
111   for (Int_t i=0; i<kMaxLayer; i++) {
112     Int_t jmax = j + fLayers[i].GetNladders()*fLayers[i].GetNdetectors();
113     for (; j<jmax; j++) {           
114       if (!cTree->GetEvent(j)) continue;
115       Int_t ncl=clusters->GetEntriesFast();
116       while (ncl--) {
117         AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
118
119 #ifdef DEBUG
120 if (c->GetLabel(2)!=LAB)
121 if (c->GetLabel(1)!=LAB)
122 if (c->GetLabel(0)!=LAB) continue;
123 cout<<lay-1<<' '<<lad-1<<' '<<det-1<<' '<<c->GetY()<<' '<<c->GetZ()<<endl;
124 #endif
125
126         fLayers[i].InsertCluster(new AliITSclusterV2(*c));
127       }
128       clusters->Delete();
129     }
130   }
131   delete cTree; //Thanks to Mariana Bondila
132 }
133
134 void AliITStrackerV2::UnloadClusters() {
135   //--------------------------------------------------------------------
136   //This function unloads ITS clusters
137   //--------------------------------------------------------------------
138   for (Int_t i=0; i<kMaxLayer; i++) fLayers[i].ResetClusters();
139 }
140
141 //#ifdef DEBUG
142 static Int_t lbl;
143 //#endif
144
145 Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) {
146   //--------------------------------------------------------------------
147   //This functions reconstructs ITS tracks
148   //--------------------------------------------------------------------
149   TFile *in=(TFile*)inp;
150   TDirectory *savedir=gDirectory; 
151
152   LoadClusters();
153
154   if (!in->IsOpen()) {
155      cerr<<"AliITStrackerV2::Clusters2Tracks(): ";
156      cerr<<"file with TPC tracks is not open !\n";
157      return 1;
158   }
159
160   if (!out->IsOpen()) {
161      cerr<<"AliITStrackerV2::Clusters2Tracks(): ";
162      cerr<<"file for ITS tracks is not open !\n";
163      return 2;
164   }
165
166   in->cd();
167  
168   Char_t tname[100];
169   Int_t nentr=0; TObjArray itsTracks(15000);
170
171   {/* Read TPC tracks */ 
172     sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
173     TTree *tpcTree=(TTree*)in->Get(tname);
174     if (!tpcTree) {
175        cerr<<"AliITStrackerV2::Clusters2Tracks(): "
176              "can't get a tree with TPC tracks !\n";
177        return 3;
178     }
179     AliTPCtrack *itrack=new AliTPCtrack; 
180     tpcTree->SetBranchAddress("tracks",&itrack);
181     nentr=(Int_t)tpcTree->GetEntries();
182     for (Int_t i=0; i<nentr; i++) {
183        tpcTree->GetEvent(i);
184        AliITStrackV2 *t=0;
185        try {
186            t=new AliITStrackV2(*itrack);
187        } catch (const Char_t *msg) {
188            cerr<<msg<<endl;
189            delete t;
190            continue;
191        }
192        if (TMath::Abs(t->GetD())>4) continue;
193
194        t->PropagateTo(80.,0.0053,30);
195        if (TMath::Abs(t->GetY())>12.8) t->CorrectForMaterial(0.03);
196        if (TMath::Abs(t->GetZ())<0.2) t->CorrectForMaterial(0.40);
197        t->PropagateTo(61.,0.0053,30);
198        //Double_t xk=52.,x,y,z; t->GetGlobalXYZat(xk,x,y,z);
199        //if (TMath::Abs(y)<7.77) t->PropagateTo(xk,0.19,24.); 
200        t->PropagateTo(50.,0.001);
201
202        itsTracks.AddLast(t);
203     }
204     delete tpcTree; //Thanks to Mariana Bondila
205     delete itrack;
206   }
207   itsTracks.Sort();
208
209   out->cd();
210
211   sprintf(tname,"TreeT_ITS_%d",GetEventNumber());
212   TTree itsTree(tname,"Tree with ITS tracks");
213   AliITStrackV2 *otrack=&fBestTrack;
214   itsTree.Branch("tracks","AliITStrackV2",&otrack,32000,0);
215
216   for (fPass=0; fPass<2; fPass++) {
217      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
218      for (Int_t i=0; i<nentr; i++) {
219        if (i%10==0) cerr<<nentr-i<<" \r";
220        AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i);
221        if (t==0) continue;           //this track has been already tracked
222        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
223
224 lbl=tpcLabel;
225 #ifdef DEBUG
226 lbl=tpcLabel;
227 if (TMath::Abs(tpcLabel)!=LAB) continue;
228 cout<<tpcLabel<<" *****************\n";
229 #endif
230
231        ResetTrackToFollow(*t);
232        ResetBestTrack();
233
234        for (FollowProlongation(); fI<kMaxLayer; fI++) {
235           while (TakeNextProlongation()) FollowProlongation();
236        }
237
238 #ifdef DEBUG
239 cout<<fBestTrack.GetNumberOfClusters()<<" number of clusters\n\n";
240 #endif
241
242       if (fBestTrack.GetNumberOfClusters() < kMaxLayer-kLayersToSkip) continue;
243
244         if (fConstraint[fPass]) {
245            Int_t index[kMaxLayer];
246            Int_t k;
247            for (k=0; k<kMaxLayer; k++) index[k]=-1;
248            Int_t nc=fBestTrack.GetNumberOfClusters();
249            for (k=0; k<nc; k++) { 
250               Int_t idx=fBestTrack.GetClusterIndex(k),nl=(idx&0xf0000000)>>28;
251               index[nl]=idx; 
252            }
253            fBestTrack.~AliITStrackV2(); new(&fBestTrack) AliITStrackV2(*t);
254            if (!RefitAt(3.7, &fBestTrack, index)) continue;
255         }
256
257         fBestTrack.SetLabel(tpcLabel);
258         fBestTrack.CookdEdx();
259         CookLabel(&fBestTrack,0.); //For comparison only
260         itsTree.Fill();
261         UseClusters(&fBestTrack);
262         delete itsTracks.RemoveAt(i);
263
264      }
265   }
266
267   itsTree.Write();
268
269   itsTracks.Delete();
270
271   UnloadClusters();
272
273   savedir->cd();
274   cerr<<"Number of TPC tracks: "<<nentr<<endl;
275   cerr<<"Number of prolonged tracks: "<<itsTree.GetEntries()<<endl;
276
277   return 0;
278 }
279
280 Int_t AliITStrackerV2::PropagateBack(const TFile *inp, TFile *out) {
281   //--------------------------------------------------------------------
282   //This functions propagates reconstructed ITS tracks back
283   //--------------------------------------------------------------------
284   TFile *in=(TFile*)inp;
285   TDirectory *savedir=gDirectory; 
286
287   LoadClusters();
288
289   if (!in->IsOpen()) {
290      cerr<<"AliITStrackerV2::PropagateBack(): ";
291      cerr<<"file with ITS tracks is not open !\n";
292      return 1;
293   }
294
295   if (!out->IsOpen()) {
296      cerr<<"AliITStrackerV2::PropagateBack(): ";
297      cerr<<"file for back propagated ITS tracks is not open !\n";
298      return 2;
299   }
300
301   in->cd();
302
303   Char_t tname[100];
304   sprintf(tname,"TreeT_ITS_%d",GetEventNumber());
305   TTree *itsTree=(TTree*)in->Get(tname);
306   if (!itsTree) {
307      cerr<<"AliITStrackerV2::PropagateBack() ";
308      cerr<<"can't get a tree with ITS tracks !\n";
309      return 3;
310   }
311   AliITStrackV2 *itrack=new AliITStrackV2; 
312   itsTree->SetBranchAddress("tracks",&itrack);
313
314   out->cd();
315
316   sprintf(tname,"TreeT_ITSb_%d",GetEventNumber());
317   TTree backTree(tname,"Tree with back propagated ITS tracks");
318   AliTPCtrack *otrack=0;
319   backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
320
321   Int_t ntrk=0;
322
323   Int_t nentr=(Int_t)itsTree->GetEntries();
324   for (Int_t i=0; i<nentr; i++) {
325     itsTree->GetEvent(i);
326     ResetTrackToFollow(*itrack);
327     fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters();
328     Int_t itsLabel=fTrackToFollow.GetLabel(); //save the ITS track label
329
330 #ifdef DEBUG
331 if (TMath::Abs(itsLabel)!=LAB) continue;
332 cout<<itsLabel<<" *****************\n";
333 #endif
334
335     try {
336        Int_t nc=itrack->GetNumberOfClusters();
337 #ifdef DEBUG
338 for (Int_t k=0; k<nc; k++) {
339     Int_t index=itrack->GetClusterIndex(k);
340     AliITSclusterV2 *c=(AliITSclusterV2*)GetCluster(index);
341     cout<<c->GetLabel(0)<<' '<<c->GetLabel(1)<<' '<<c->GetLabel(2)<<endl;
342 }
343 #endif       
344        Int_t idx=0, l=0; 
345        const  AliITSclusterV2 *c=0; 
346        if (nc--) {
347           idx=itrack->GetClusterIndex(nc); l=(idx&0xf0000000)>>28;
348           c=(AliITSclusterV2*)GetCluster(idx);
349        }
350        for (fI=0; fI<kMaxLayer; fI++) {
351          AliITSlayer &layer=fLayers[fI];
352          Double_t r=layer.GetR();
353          if (fI==2 || fI==4) {             
354             Double_t rs=0.5*(fLayers[fI-1].GetR() + r);
355             Double_t d=0.0034, x0=38.6;
356             if (fI==2) {rs=9.; d=0.0097; x0=42.;}
357             if (!fTrackToFollow.PropagateTo(rs,-d,x0)) throw "";
358          }
359
360          Double_t x,y,z;
361          if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) 
362             throw "AliITStrackerV2::PropagateBack: failed to estimate track !";
363          Double_t phi=TMath::ATan2(y,x);
364          Int_t idet=layer.FindDetectorIndex(phi,z);
365          if (idet<0) 
366          throw "AliITStrackerV2::PropagateBack: failed to find a detector !\n";
367          const AliITSdetector &det=layer.GetDetector(idet);
368          r=det.GetR(); phi=det.GetPhi();
369          if (!fTrackToFollow.Propagate(phi,r)) throw "";
370          fTrackToFollow.SetDetectorIndex(idet);
371
372          const AliITSclusterV2 *cl=0;
373          Int_t index=0;
374          Double_t maxchi2=kMaxChi2;
375
376          if (l==fI) {
377            idet=c->GetDetectorIndex();
378            if (idet != fTrackToFollow.GetDetectorIndex()) {
379              const AliITSdetector &det=layer.GetDetector(idet);
380              r=det.GetR(); phi=det.GetPhi();
381              if (!fTrackToFollow.Propagate(phi,r)) throw "";
382              fTrackToFollow.SetDetectorIndex(idet);
383            }
384            Double_t chi2=fTrackToFollow.GetPredictedChi2(c);
385            if (chi2<kMaxChi2) {
386               cl=c; maxchi2=chi2; index=idx;
387            }
388            if (nc--) {
389               idx=itrack->GetClusterIndex(nc); l=(idx&0xf0000000)>>28;
390               c=(AliITSclusterV2*)GetCluster(idx);
391            } 
392          }
393
394          if (fTrackToFollow.GetNumberOfClusters()>2) {
395            Double_t dz=3*TMath::Sqrt(fTrackToFollow.GetSigmaZ2()+kSigmaZ2[fI]);
396            Double_t dy=3*TMath::Sqrt(fTrackToFollow.GetSigmaY2()+kSigmaY2[fI]);
397            Double_t zmin=fTrackToFollow.GetZ() - dz;
398            Double_t zmax=fTrackToFollow.GetZ() + dz;
399            Double_t ymin=fTrackToFollow.GetY() + phi*r - dy;
400            Double_t ymax=fTrackToFollow.GetY() + phi*r + dy;
401            layer.SelectClusters(zmin,zmax,ymin,ymax);
402
403            const AliITSclusterV2 *cc=0; Int_t ci;
404            while ((cc=layer.GetNextCluster(ci))!=0) {
405               idet=cc->GetDetectorIndex();
406               if (idet != fTrackToFollow.GetDetectorIndex()) continue;
407               Double_t chi2=fTrackToFollow.GetPredictedChi2(cc);
408               if (chi2<maxchi2) {
409                  cl=cc; index=(fI<<28)+ci; maxchi2=chi2;
410               }
411            }
412          }
413
414          if (cl) {
415             if (!fTrackToFollow.Update(cl,maxchi2,index)) 
416               cerr<<"AliITStrackerV2::PropagateBack: filtering failed !\n";
417          }
418          {
419          Double_t x0;
420          x=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ(),x0);
421          fTrackToFollow.CorrectForMaterial(-x,x0); 
422          }
423        }
424
425        fTrackToFollow.PropagateTo(50.,-0.001);
426        //Double_t xk=52.,x,y,z; fTrackToFollow.GetGlobalXYZat(xk,x,y,z);
427        //if (TMath::Abs(y)<7.77) fTrackToFollow.PropagateTo(xk,-0.19,24.); 
428        fTrackToFollow.PropagateTo(61,-0.0053,30);
429        fTrackToFollow.PropagateTo(80.,-0.0053,30);
430
431        fTrackToFollow.SetLabel(itsLabel);
432        otrack=new AliTPCtrack(fTrackToFollow,fTrackToFollow.GetAlpha()); 
433        backTree.Fill(); delete otrack;
434        UseClusters(&fTrackToFollow);
435        cerr << ++ntrk << "                \r";
436     }
437     catch (const Char_t *msg) {
438        cerr<<msg<<endl;
439     }
440   }
441
442   backTree.Write();
443
444   cerr<<"Number of ITS tracks: "<<nentr<<endl;
445   cerr<<"Number of back propagated ITS tracks: "<<ntrk<<endl;
446
447   delete itrack;
448
449   delete itsTree; //Thanks to Mariana Bondila
450
451   UnloadClusters();
452
453   savedir->cd();
454
455   return 0;
456 }
457
458 AliCluster *AliITStrackerV2::GetCluster(Int_t index) const {
459   //--------------------------------------------------------------------
460   //       Return pointer to a given cluster
461   //--------------------------------------------------------------------
462   Int_t l=(index & 0xf0000000) >> 28;
463   Int_t c=(index & 0x0fffffff) >> 00;
464   return fLayers[l].GetCluster(c);
465 }
466
467
468 void AliITStrackerV2::FollowProlongation() {
469   //--------------------------------------------------------------------
470   //This function finds a track prolongation 
471   //--------------------------------------------------------------------
472   Int_t tryAgain=kLayersToSkip;
473
474   while (fI) {
475     Int_t i=fI-1;
476 #ifdef DEBUG
477 cout<<i<<' ';
478 #endif
479     AliITSlayer &layer=fLayers[i];
480     AliITStrackV2 &track=fTracks[i];
481
482     Double_t r=layer.GetR();
483
484     if (i==3 || i==1) {
485        Double_t rs=0.5*(fLayers[i+1].GetR() + r);
486        Double_t d=0.0034, x0=38.6;
487        if (i==1) {rs=9.; d=0.0097; x0=42;}
488        if (!fTrackToFollow.PropagateTo(rs,d,x0)) {
489          //cerr<<"AliITStrackerV2::FollowProlongation: "
490          //"propagation failed !\n";
491          break;
492        }
493     }
494
495     //find intersection
496     Double_t x,y,z;  
497     if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) {
498       //cerr<<"AliITStrackerV2::FollowProlongation: "
499       //"failed to estimate track !\n";
500       break;
501     }
502     Double_t phi=TMath::ATan2(y,x);
503     Int_t idet=layer.FindDetectorIndex(phi,z);
504     if (idet<0) {
505       //cerr<<"AliITStrackerV2::FollowProlongation: "
506       //"failed to find a detector !\n";
507       break;
508     }
509
510     //propagate to the intersection
511     const AliITSdetector &det=layer.GetDetector(idet);
512     phi=det.GetPhi();
513     if (!fTrackToFollow.Propagate(phi,det.GetR())) {
514       //cerr<<"AliITStrackerV2::FollowProlongation: propagation failed !\n";
515       break;
516     }
517     fTrackToFollow.SetDetectorIndex(idet);
518
519     //Select possible prolongations and store the current track estimation
520     track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow);
521     Double_t dz=3*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]);
522     if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl());
523     if (dz > kMaxRoad) {
524       //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Z !\n";
525       break;
526     }
527
528     if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+dz) break;
529
530     Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]);
531     if (dy < 0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp());
532     if (dy > kMaxRoad) {
533       //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Y !\n";
534       break;
535     }
536
537     Double_t zmin=track.GetZ() - dz; 
538     Double_t zmax=track.GetZ() + dz;
539     Double_t ymin=track.GetY() + r*phi - dy;
540     Double_t ymax=track.GetY() + r*phi + dy;
541     layer.SelectClusters(zmin,zmax,ymin,ymax); 
542     fI--;
543
544     //take another prolongation
545     if (!TakeNextProlongation()) if (!tryAgain--) break;
546     tryAgain=kLayersToSkip;
547
548   } 
549
550   //deal with the best track
551   Int_t ncl=fTrackToFollow.GetNumberOfClusters();
552   Int_t nclb=fBestTrack.GetNumberOfClusters();
553   if (ncl)
554   if (ncl >= nclb) {
555      Double_t chi2=fTrackToFollow.GetChi2();
556      if (chi2/ncl < kChi2PerCluster) {        
557         if (ncl > nclb || chi2 < fBestTrack.GetChi2()) {
558            ResetBestTrack();
559         }
560      }
561   }
562
563 }
564
565 Int_t AliITStrackerV2::TakeNextProlongation() {
566   //--------------------------------------------------------------------
567   // This function takes another track prolongation 
568   //
569   //  dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch 
570   //--------------------------------------------------------------------
571   AliITSlayer &layer=fLayers[fI];
572   ResetTrackToFollow(fTracks[fI]);
573
574   Double_t dz=4*TMath::Sqrt(fTrackToFollow.GetSigmaZ2() + kSigmaZ2[fI]);
575   Double_t dy=4*TMath::Sqrt(fTrackToFollow.GetSigmaY2() + kSigmaY2[fI]);
576
577   const AliITSclusterV2 *c=0; Int_t ci=-1;
578   Double_t chi2=12345.;
579   while ((c=layer.GetNextCluster(ci))!=0) {
580     //if (c->GetLabel(0)!=TMath::Abs(lbl)) continue; 
581     Int_t idet=c->GetDetectorIndex();
582
583     if (fTrackToFollow.GetDetectorIndex()!=idet) {
584        const AliITSdetector &det=layer.GetDetector(idet);
585        ResetTrackToFollow(fTracks[fI]);
586        if (!fTrackToFollow.Propagate(det.GetPhi(),det.GetR())) {
587          //cerr<<"AliITStrackerV2::TakeNextProlongation: "
588          //"propagation failed !\n";
589          continue;
590        }
591        fTrackToFollow.SetDetectorIndex(idet);
592        if (TMath::Abs(fTrackToFollow.GetZ()-GetZ())>layer.GetR()+dz) continue;
593
594 #ifdef DEBUG
595 cout<<fI<<" change detector !\n";
596 #endif
597
598     }
599
600     if (TMath::Abs(fTrackToFollow.GetZ() - c->GetZ()) > dz) continue;
601     if (TMath::Abs(fTrackToFollow.GetY() - c->GetY()) > dy) continue;
602
603     chi2=fTrackToFollow.GetPredictedChi2(c); if (chi2<kMaxChi2) break;
604   }
605
606 #ifdef DEBUG
607 cout<<fI<<" chi2="<<chi2<<' '
608     <<fTrackToFollow.GetY()<<' '<<fTrackToFollow.GetZ()<<' '
609     <<dy<<' '<<dz<<endl;
610 {
611 Double_t phi=TMath::ATan(fTrackToFollow.GetY()/fTrackToFollow.GetX())
612             +fTrackToFollow.GetAlpha(); phi=phi*180./3.141593;
613 cout<<phi<<endl;
614 }
615 #endif
616
617   if (chi2>=kMaxChi2) return 0;
618   if (!c) return 0;
619
620   if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) {
621      //cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
622      return 0;
623   }
624
625   if (fTrackToFollow.GetNumberOfClusters()>1)
626   if (TMath::Abs(fTrackToFollow.GetD())>4) return 0;
627
628   fTrackToFollow.
629     SetSampledEdx(c->GetQ(),fTrackToFollow.GetNumberOfClusters()-1); //b.b.
630
631   {
632  Double_t x0;
633  Double_t d=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ(),x0);
634    fTrackToFollow.CorrectForMaterial(d,x0);
635   }
636
637   if (fConstraint[fPass]) {
638     Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
639     fTrackToFollow.Improve(d,GetY(),GetZ());
640   }
641
642 #ifdef DEBUG
643 cout<<"accepted lab="<<c->GetLabel(0)<<' '
644     <<fTrackToFollow.GetNumberOfClusters()<<' '
645     <<fTrackToFollow.GetY()<<' '<<fTrackToFollow.GetZ()<<' '
646     <<fTrackToFollow.Get1Pt()<<endl<<endl;
647 #endif
648
649   return 1;
650 }
651
652
653 AliITStrackerV2::AliITSlayer::AliITSlayer() {
654   //--------------------------------------------------------------------
655   //default AliITSlayer constructor
656   //--------------------------------------------------------------------
657   fN=0;
658   fDetectors=0;
659 }
660
661 AliITStrackerV2::AliITSlayer::
662 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd) {
663   //--------------------------------------------------------------------
664   //main AliITSlayer constructor
665   //--------------------------------------------------------------------
666   fR=r; fPhiOffset=p; fZOffset=z;
667   fNladders=nl; fNdetectors=nd;
668   fDetectors=new AliITSdetector[fNladders*fNdetectors];
669
670   fN=0;
671   fI=0;
672 }
673
674 AliITStrackerV2::AliITSlayer::~AliITSlayer() {
675   //--------------------------------------------------------------------
676   // AliITSlayer destructor
677   //--------------------------------------------------------------------
678   delete[] fDetectors;
679   for (Int_t i=0; i<fN; i++) delete fClusters[i];
680 }
681
682 void AliITStrackerV2::AliITSlayer::ResetClusters() {
683   //--------------------------------------------------------------------
684   // This function removes loaded clusters
685   //--------------------------------------------------------------------
686   for (Int_t i=0; i<fN; i++) delete fClusters[i];
687   fN=0;
688   fI=0;
689 }
690
691 Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSclusterV2 *c) {
692   //--------------------------------------------------------------------
693   //This function adds a cluster to this layer
694   //--------------------------------------------------------------------
695   if (fN==kMaxClusterPerLayer) {
696      cerr<<"AliITStrackerV2::AliITSlayer::InsertCluster(): "
697            "Too many clusters !\n";
698      return 1;
699   }
700
701   if (fN==0) {fClusters[fN++]=c; return 0;}
702   Int_t i=FindClusterIndex(c->GetZ());
703   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliITSclusterV2*));
704   fClusters[i]=c; fN++;
705
706   return 0;
707 }
708
709 Int_t AliITStrackerV2::AliITSlayer::FindClusterIndex(Double_t z) const {
710   //--------------------------------------------------------------------
711   // This function returns the index of the nearest cluster 
712   //--------------------------------------------------------------------
713   if (fN==0) return 0;
714   if (z <= fClusters[0]->GetZ()) return 0;
715   if (z > fClusters[fN-1]->GetZ()) return fN;
716   Int_t b=0, e=fN-1, m=(b+e)/2;
717   for (; b<e; m=(b+e)/2) {
718     if (z > fClusters[m]->GetZ()) b=m+1;
719     else e=m; 
720   }
721   return m;
722 }
723
724 void AliITStrackerV2::AliITSlayer::
725 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
726   //--------------------------------------------------------------------
727   // This function sets the "window"
728   //--------------------------------------------------------------------
729   fI=FindClusterIndex(zmin); fZmax=zmax;
730   Double_t circle=2*TMath::Pi()*fR;
731   if (ymax>circle) { ymax-=circle; ymin-=circle; }
732   fYmin=ymin; fYmax=ymax;
733 }
734
735 const AliITSclusterV2 *AliITStrackerV2::AliITSlayer::GetNextCluster(Int_t &ci){
736   //--------------------------------------------------------------------
737   // This function returns clusters within the "window" 
738   //--------------------------------------------------------------------
739   const AliITSclusterV2 *cluster=0;
740   for (Int_t i=fI; i<fN; i++) {
741     const AliITSclusterV2 *c=fClusters[i];
742     if (c->GetZ() > fZmax) break;
743     if (c->IsUsed()) continue;
744     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());    
745     Double_t y=fR*det.GetPhi() + c->GetY();
746
747     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
748     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
749
750     if (y<fYmin) continue;
751     if (y>fYmax) continue;
752     cluster=c; ci=i;
753     fI=i+1;
754     break; 
755   }
756
757   return cluster;
758 }
759
760 Int_t AliITStrackerV2::AliITSlayer::
761 FindDetectorIndex(Double_t phi, Double_t z) const {
762   //--------------------------------------------------------------------
763   //This function finds the detector crossed by the track
764   //--------------------------------------------------------------------
765   Double_t dphi=phi-fPhiOffset;
766   if      (dphi <  0) dphi += 2*TMath::Pi();
767   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
768   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
769   if (np>=fNladders) np-=fNladders;
770   if (np<0)          np+=fNladders;
771
772   Double_t dz=fZOffset-z;
773   Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
774   if (nz>=fNdetectors) return -1;
775   if (nz<0)            return -1;
776
777 #ifdef DEBUG
778 cout<<np<<' '<<nz<<endl;
779 #endif
780
781   return np*fNdetectors + nz;
782 }
783
784 Double_t 
785 AliITStrackerV2::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
786 const {
787   //--------------------------------------------------------------------
788   //This function returns the layer thickness at this point (units X0)
789   //--------------------------------------------------------------------
790   Double_t d=0.0085;
791   x0=21.82;
792
793   if (43<fR&&fR<45) { //SSD2
794      Double_t dd=0.0034;
795      d=dd;
796      if (TMath::Abs(y-0.00)>3.40) d+=dd;
797      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
798      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
799      for (Int_t i=0; i<12; i++) {
800        if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
801           if (TMath::Abs(y-0.00)>3.40) d+=dd;
802           d+=0.0034; 
803           break;
804        }
805        if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
806           if (TMath::Abs(y-0.00)>3.40) d+=dd;
807           d+=0.0034; 
808           break;
809        }         
810        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
811        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
812      }
813   } else 
814   if (37<fR&&fR<41) { //SSD1
815      Double_t dd=0.0034;
816      d=dd;
817      if (TMath::Abs(y-0.00)>3.40) d+=dd;
818      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
819      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
820      for (Int_t i=0; i<11; i++) {
821        if (TMath::Abs(z-3.9*i)<0.15) {
822           if (TMath::Abs(y-0.00)>3.40) d+=dd;
823           d+=dd; 
824           break;
825        }
826        if (TMath::Abs(z+3.9*i)<0.15) {
827           if (TMath::Abs(y-0.00)>3.40) d+=dd;
828           d+=dd; 
829           break;
830        }         
831        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
832        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}         
833      }
834   } else
835   if (13<fR&&fR<26) { //SDD
836      Double_t dd=0.0033;
837      d=dd;
838      if (TMath::Abs(y-0.00)>3.30) d+=dd;
839
840      if (TMath::Abs(y-1.80)<0.55) {
841         d+=0.016;
842         for (Int_t j=0; j<20; j++) {
843           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
844           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
845         } 
846      }
847      if (TMath::Abs(y+1.80)<0.55) {
848         d+=0.016;
849         for (Int_t j=0; j<20; j++) {
850           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
851           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
852         } 
853      }
854
855      for (Int_t i=0; i<4; i++) {
856        if (TMath::Abs(z-7.3*i)<0.60) {
857           d+=dd;
858           if (TMath::Abs(y-0.00)>3.30) d+=dd; 
859           break;
860        }
861        if (TMath::Abs(z+7.3*i)<0.60) {
862           d+=dd; 
863           if (TMath::Abs(y-0.00)>3.30) d+=dd; 
864           break;
865        }
866      }
867   } else
868   if (6<fR&&fR<8) {   //SPD2
869      Double_t dd=0.0063; x0=21.5;
870      d=dd;
871      if (TMath::Abs(y-3.08)>0.5) d+=dd;
872      //if (TMath::Abs(y-3.08)>0.45) d+=dd;
873      if (TMath::Abs(y-3.03)<0.10) {d+=0.014;}
874   } else
875   if (3<fR&&fR<5) {   //SPD1
876      Double_t dd=0.0063; x0=21.5;
877      d=dd;
878      if (TMath::Abs(y+0.21)>0.6) d+=dd;
879      //if (TMath::Abs(y+0.21)>0.45) d+=dd;
880      if (TMath::Abs(y+0.10)<0.10) {d+=0.014;}
881   }
882
883 #ifdef DEBUG
884   cout<<"d="<<d<<endl;
885 #endif
886
887   return d;
888 }
889
890 Double_t AliITStrackerV2::GetEffectiveThickness(Double_t y,Double_t z) const
891 {
892   //--------------------------------------------------------------------
893   //Returns the thickness between the current layer and the vertex (units X0)
894   //--------------------------------------------------------------------
895   Double_t d=0.0028*3*3; //beam pipe
896   Double_t x0=0;
897
898   Double_t xn=fLayers[fI].GetR();
899   for (Int_t i=0; i<fI; i++) {
900     Double_t xi=fLayers[i].GetR();
901     d+=fLayers[i].GetThickness(y,z,x0)*xi*xi;
902   }
903
904   if (fI>1) {
905     Double_t xi=9.;
906     d+=0.0097*xi*xi;
907   }
908
909   if (fI>3) {
910     Double_t xi=0.5*(fLayers[3].GetR()+fLayers[4].GetR());
911     d+=0.0034*xi*xi;
912   }
913
914   return d/(xn*xn);
915 }
916
917 Int_t AliITStrackerV2::AliITSlayer::InRoad() const {
918   //--------------------------------------------------------------------
919   // This function returns number of clusters within the "window" 
920   //--------------------------------------------------------------------
921   Int_t ncl=0;
922   for (Int_t i=fI; i<fN; i++) {
923     const AliITSclusterV2 *c=fClusters[i];
924     if (c->GetZ() > fZmax) break;
925     if (c->IsUsed()) continue;
926     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());    
927     Double_t y=fR*det.GetPhi() + c->GetY();
928
929     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
930     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
931
932     if (y<fYmin) continue;
933     if (y>fYmax) continue;
934     ncl++;
935   }
936   return ncl;
937 }
938
939 Bool_t AliITStrackerV2::RefitAt(Double_t x, AliITStrackV2 *t, Int_t *index) {
940   //--------------------------------------------------------------------
941   // This function refits a track at a given position
942   //--------------------------------------------------------------------
943   Int_t from, to, step;
944   if (x > t->GetX()) {
945       from=0; to=kMaxLayer;
946       step=+1;
947   } else {
948       from=kMaxLayer-1; to=-1;
949       step=-1;
950   }
951
952   for (Int_t i=from; i != to; i += step) {
953      AliITSlayer &layer=fLayers[i];
954      Double_t r=layer.GetR();
955  
956      {
957      Double_t hI=i-0.5*step; 
958      if (hI==1.5 || hI==3.5) {             
959         Double_t rs=0.5*(fLayers[i-step].GetR() + r);
960         Double_t d=0.0034, x0=38.6; 
961         if (hI==1.5) {rs=9.; d=0.0097; x0=42;}
962         if (!t->PropagateTo(rs,d,x0)) {
963           return kFALSE;
964         }
965      }
966      }
967
968      Double_t x,y,z;
969      if (!t->GetGlobalXYZat(r,x,y,z)) { 
970        return kFALSE;
971      }
972      Double_t phi=TMath::ATan2(y,x);
973      Int_t idet=layer.FindDetectorIndex(phi,z);
974      if (idet<0) { 
975        return kFALSE;
976      }
977      const AliITSdetector &det=layer.GetDetector(idet);
978      phi=det.GetPhi();
979      if (!t->Propagate(phi,det.GetR())) {
980        return kFALSE;
981      }
982      t->SetDetectorIndex(idet);
983
984      const AliITSclusterV2 *cl=0;
985      Double_t maxchi2=kMaxChi2;
986
987      Int_t idx=index[i];
988      if (idx>0) {
989         const AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(idx); 
990         if (idet != c->GetDetectorIndex()) {
991            idet=c->GetDetectorIndex();
992            const AliITSdetector &det=layer.GetDetector(idet);
993            if (!t->Propagate(det.GetPhi(),det.GetR())) {
994              return kFALSE;
995            }
996            t->SetDetectorIndex(idet);
997         }
998         Double_t chi2=t->GetPredictedChi2(c);
999         if (chi2<maxchi2) { cl=c; maxchi2=chi2; }
1000         else return kFALSE;
1001      }
1002
1003      /*
1004      if (cl==0)
1005      if (t->GetNumberOfClusters()>2) {
1006         Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]);
1007         Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
1008         Double_t zmin=t->GetZ() - dz;
1009         Double_t zmax=t->GetZ() + dz;
1010         Double_t ymin=t->GetY() + phi*r - dy;
1011         Double_t ymax=t->GetY() + phi*r + dy;
1012         layer.SelectClusters(zmin,zmax,ymin,ymax);
1013
1014         const AliITSclusterV2 *c=0; Int_t ci=-1;
1015         while ((c=layer.GetNextCluster(ci))!=0) {
1016            if (idet != c->GetDetectorIndex()) continue;
1017            Double_t chi2=t->GetPredictedChi2(c);
1018            if (chi2<maxchi2) { cl=c; maxchi2=chi2; idx=ci; }
1019         }
1020      }
1021      */
1022
1023      if (cl) {
1024        if (!t->Update(cl,maxchi2,idx)) {
1025           return kFALSE;
1026        }
1027      }
1028
1029      {
1030      Double_t x0;
1031      Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1032      t->CorrectForMaterial(-step*d,x0);
1033      }
1034
1035   }
1036
1037   if (!t->PropagateTo(x,0.,0.)) return kFALSE;
1038   return kTRUE;
1039 }
1040
1041
1042