First commit of ITS tracking version 2 (Yu.Belikov)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 13 Mar 2001 14:24:06 +0000 (14:24 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 13 Mar 2001 14:24:06 +0000 (14:24 +0000)
13 files changed:
ITS/AliITSComparisonV2.C [new file with mode: 0644]
ITS/AliITSFindClustersV2.C [new file with mode: 0644]
ITS/AliITSFindTracksV2.C [new file with mode: 0644]
ITS/AliITSclusterV2.cxx [new file with mode: 0644]
ITS/AliITSclusterV2.h [new file with mode: 0644]
ITS/AliITSrecoV2.h [new file with mode: 0644]
ITS/AliITStestV2.C [new file with mode: 0644]
ITS/AliITStrackV2.cxx [new file with mode: 0644]
ITS/AliITStrackV2.h [new file with mode: 0644]
ITS/AliITStrackerV2.cxx [new file with mode: 0644]
ITS/AliITStrackerV2.h [new file with mode: 0644]
ITS/ITSLinkDef.h
ITS/Makefile

diff --git a/ITS/AliITSComparisonV2.C b/ITS/AliITSComparisonV2.C
new file mode 100644 (file)
index 0000000..ee287dc
--- /dev/null
@@ -0,0 +1,326 @@
+#ifndef __CINT__
+  #include <iostream.h>
+  #include <fstream.h>
+
+  #include "AliRun.h"
+  #include "AliITS.h"
+  #include "AliITSgeom.h"
+  #include "AliITStrackerV2.h"
+  #include "AliITStrackV2.h"
+  #include "AliITSclusterV2.h"
+
+  #include "TFile.h"
+  #include "TTree.h"
+  #include "TH1.h"
+  #include "TObjArray.h"
+  #include "TStyle.h"
+  #include "TCanvas.h"
+  #include "TLine.h"
+  #include "TText.h"
+  #include "TParticle.h"
+#endif
+
+struct GoodTrack {
+  Int_t lab;
+  Int_t code;
+  Float_t px,py,pz;
+  Float_t x,y,z;
+};
+Int_t good_tracks(GoodTrack *gt, Int_t max);
+
+Int_t AliITSComparisonV2() {
+   cerr<<"Doing comparison...\n";
+
+   TFile *cf=TFile::Open("AliITSclustersV2.root");
+   if (!cf->IsOpen()) {cerr<<"Can't open AliITSclustersV2.root !\n"; return 1;}
+   AliITSgeom *geom=(AliITSgeom*)cf->Get("AliITSgeom");
+   if (!geom) { cerr<<"Can't get the ITS geometry !\n"; return 2; }
+   AliITStrackerV2 tracker(geom);   
+
+// Load tracks
+   TFile *tf=TFile::Open("AliITStracksV2.root");
+   if (!tf->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 3;}
+   TObjArray tarray(2000);
+   TTree *tracktree=(TTree*)tf->Get("TreeT");
+   TBranch *tbranch=tracktree->GetBranch("tracks");
+   Int_t nentr=(Int_t)tracktree->GetEntries(),i;
+   for (i=0; i<nentr; i++) {
+       AliITStrackV2 *iotrack=new AliITStrackV2;
+       tbranch->SetAddress(&iotrack);
+       tracktree->GetEvent(i);
+
+       Int_t tpcLabel=iotrack->GetLabel();
+       tracker.CookLabel(iotrack,0.);
+       Int_t itsLabel=iotrack->GetLabel();
+       if (itsLabel != tpcLabel) iotrack->SetLabel(-TMath::Abs(itsLabel));
+       if (tpcLabel < 0)         iotrack->SetLabel(-TMath::Abs(itsLabel));
+       /*
+       if (itsLabel==1234) {
+         Int_t nc=iotrack->GetNumberOfClusters();
+         for (Int_t k=0; k<nc; k++) {
+           Int_t index=iotrack->GetClusterIndex(k);
+           AliITSclusterV2 *c=tracker.GetCluster(index);
+           cout<<c->GetLabel(0)<<' '<<c->GetLabel(1)<<' '<<c->GetLabel(2)<<endl;
+         }
+       }
+       */
+       tarray.AddLast(iotrack);
+   }   
+   tf->Close();
+   cf->Close();
+
+/////////////////////////////////////////////////////////////////////////
+   GoodTrack gt[15000];
+   Int_t ngood=0;
+   ifstream in("good_tracks_its");
+   if (in) {
+      cerr<<"Reading good tracks...\n";
+      while (in>>gt[ngood].lab>>gt[ngood].code>>
+                 gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>>
+                 gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
+         ngood++;
+         cerr<<ngood<<'\r';
+         if (ngood==15000) {
+            cerr<<"Too many good tracks !\n";
+            break;
+         }
+      }
+      if (!in.eof()) cerr<<"Read error (good_tracks_its) !\n";
+   } else {
+      cerr<<"Marking good tracks (this will take a while)...\n";
+      ngood=good_tracks(gt,15000);
+      ofstream out("good_tracks_its");
+      if (out) {
+         for (Int_t ngd=0; ngd<ngood; ngd++)            
+           out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '<<
+                 gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '<<
+                 gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
+      } else cerr<<"Can not open file (good_tracks_its) !\n";
+      out.close();
+   }
+   cerr<<"Number of good tracks : "<<ngood<<endl;
+
+   TH1F *hp=new TH1F("hp","PHI resolution",50,-20.,20.); hp->SetFillColor(4);
+   TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-20,20);hl->SetFillColor(4);
+   TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.); 
+   hpt->SetFillColor(2); 
+   TH1F *hmpt=new TH1F("hmpt","Transverse impact parameter",30,-300,300); 
+   hmpt->SetFillColor(6);
+   TH1F *hz=new TH1F("hz","Longitudinal impact parameter",30,-300,300); 
+   //hmpt->SetFillColor(6);
+
+   TH1F *hgood=new TH1F("hgood","Good tracks",30,0.1,6.1);    
+   TH1F *hfound=new TH1F("hfound","Found tracks",30,0.1,6.1);
+   TH1F *hfake=new TH1F("hfake","Fake tracks",30,0.1,6.1);
+   TH1F *hg=new TH1F("hg","",30,0.1,6.1); //efficiency for good tracks
+   hg->SetLineColor(4); hg->SetLineWidth(2);
+   TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,0.1,6.1);
+   hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
+
+   TH1F *hptw=new TH1F("hptw","Weghted pt",30,0.1,6.1);
+
+   while (ngood--) {
+      Int_t lab=gt[ngood].lab, tlab=-1;
+      Double_t pxg=gt[ngood].px, pyg=gt[ngood].py, pzg=gt[ngood].pz;
+      Double_t ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
+
+      hgood->Fill(ptg);
+
+      AliITStrackV2 *track=0;
+      Int_t j;
+      for (j=0; j<nentr; j++) {
+          track=(AliITStrackV2*)tarray.UncheckedAt(j);
+          tlab=track->GetLabel();
+          if (lab==TMath::Abs(tlab)) break;
+      }
+      if (j==nentr) {
+       cerr<<"Track "<<lab<<" was not found !\n";
+        continue;
+      }
+      track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+      track->PropagateToVertex();
+
+      if (lab==tlab) hfound->Fill(ptg);
+      else { hfake->Fill(ptg); cerr<<lab<<" fake\n";}
+
+      Double_t xv,par[5]; track->GetExternalParameters(xv,par);
+      Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
+      if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+      if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+      Float_t lam=TMath::ATan(par[3]); 
+      Float_t pt_1=TMath::Abs(par[4]);
+
+      Double_t phig=TMath::ATan2(pyg,pxg);
+      hp->Fill((phi - phig)*1000.);
+
+      Double_t lamg=TMath::ATan2(pzg,ptg);
+      hl->Fill((lam - lamg)*1000.);
+
+      Double_t d=10000*track->GetD();
+      hmpt->Fill(d);
+
+      hptw->Fill(ptg,TMath::Abs(d));
+
+      Double_t z=10000*track->GetZ();
+      hz->Fill(z);
+
+//if (TMath::Abs(gt[ngood].code)==11 && ptg>4.)
+      hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
+
+   }
+
+   Stat_t ng=hgood->GetEntries(); cerr<<"Good tracks "<<ng<<endl;
+   Stat_t nf=hfound->GetEntries();
+   if (ng!=0) 
+      cerr<<"Integral efficiency is about "<<nf/ng*100.<<" %\n";
+
+   gStyle->SetOptStat(111110);
+   gStyle->SetOptFit(1);
+
+   TCanvas *c1=new TCanvas("c1","",0,0,700,850);
+
+   TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
+   p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10); 
+   hp->SetFillColor(4);  hp->SetXTitle("(mrad)"); hp->Fit("gaus"); c1->cd();
+
+   TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw(); 
+   p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
+   hl->SetXTitle("(mrad)"); hl->Fit("gaus"); c1->cd();
+
+   TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
+   p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10); 
+   hpt->SetXTitle("(%)"); hpt->Fit("gaus"); c1->cd();
+
+   TPad *p4=new TPad("p4","",0.5,0,1,0.3); p4->Draw();
+   p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
+   hmpt->SetXTitle("(micron)"); hmpt->Fit("gaus"); hz->Draw("same"); c1->cd();
+   //hfound->Sumw2();
+   //hptw->Sumw2(); 
+   //hg->SetMaximum(333);
+   //hg->SetYTitle("Impact Parameter Resolution (micron)");
+   //hg->SetXTitle("Pt (GeV/c)");
+   //hg->GetXaxis()->SetRange(0,10);
+   //hg->Divide(hptw,hfound,1,1.);
+   //hg->DrawCopy(); c1->cd();
+   
+
+   TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd(); 
+   p5->SetFillColor(41); p5->SetFrameFillColor(10);
+   hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
+   hg->Divide(hfound,hgood,1,1.,"b");
+   hf->Divide(hfake,hgood,1,1.,"b");
+   hg->SetMaximum(1.4);
+   hg->SetYTitle("Tracking efficiency");
+   hg->SetXTitle("Pt (GeV/c)");
+   hg->Draw();
+
+   TLine *line1 = new TLine(0,1.0,7,1.0); line1->SetLineStyle(4);
+   line1->Draw("same");
+   TLine *line2 = new TLine(0,0.9,7,0.9); line2->SetLineStyle(4);
+   line2->Draw("same");
+
+   hf->SetFillColor(1);
+   hf->SetFillStyle(3013);
+   hf->SetLineColor(2);
+   hf->SetLineWidth(2);
+   hf->Draw("histsame");
+   TText *text = new TText(0.461176,0.248448,"Fake tracks");
+   text->SetTextSize(0.05);
+   text->Draw();
+   text = new TText(0.453919,1.11408,"Good tracks");
+   text->SetTextSize(0.05);
+   text->Draw();
+
+   return 0;
+}
+
+Int_t good_tracks(GoodTrack *gt, Int_t max) {
+   if (gAlice) {delete gAlice; gAlice=0;}
+
+   TFile *file=TFile::Open("galice.root");
+   if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
+   if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
+     cerr<<"gAlice have not been found on galice.root !\n";
+     exit(5);
+   }
+
+   Int_t np=gAlice->GetEvent(0);
+
+   Int_t *good=new Int_t[np];
+   Int_t k;
+   for (k=0; k<np; k++) good[k]=0;
+
+   AliITS *ITS=(AliITS*)gAlice->GetDetector("ITS");
+   if (!ITS) {
+      cerr<<"can't get ITS !\n"; exit(8);
+   }
+   AliITSgeom *geom=ITS->GetITSgeom();
+   if (!geom) {
+      cerr<<"cen't get ITS geometry !\n"; exit(9);
+   }
+
+   TFile *cf=TFile::Open("AliITSclustersV2.root");
+   if (!cf->IsOpen()){
+      cerr<<"Can't open AliITSclustersV2.root !\n"; exit(6);
+   }
+   TTree *cTree=(TTree*)cf->Get("cTree");
+   if (!cTree) {
+      cerr<<"Can't get cTree !\n"; exit(7);
+   }
+   TBranch *branch=cTree->GetBranch("Clusters");
+   if (!branch) {
+      cerr<<"Can't get clusters branch !\n"; exit(8);
+   }
+   TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
+   branch->SetAddress(&clusters);
+
+   Int_t entr=(Int_t)cTree->GetEntries();
+   for (k=0; k<entr; k++) {
+     if (!cTree->GetEvent(k)) continue;
+     Int_t lay,lad,det;  geom->GetModuleId(k-1,lay,lad,det);
+     if (lay<1 || lay>6) {
+       cerr<<"wrong layer !\n"; exit(10);
+     }
+     Int_t ncl=clusters->GetEntriesFast();
+     while (ncl--) {
+        AliITSclusterV2 *pnt=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
+        Int_t l0=pnt->GetLabel(0);
+        Int_t l1=pnt->GetLabel(1);
+        Int_t l2=pnt->GetLabel(2);
+        Int_t mask=1<<(lay-1);
+        if (l0>=0) good[l0]|=mask; 
+        if (l1>=0) good[l1]|=mask; 
+        if (l2>=0) good[l2]|=mask;
+     }
+   }
+   clusters->Delete(); delete clusters;
+   cf->Close();
+
+   ifstream in("good_tracks_tpc");
+   if (!in) {
+     cerr<<"can't get good_tracks_tpc !\n"; exit(11);
+   }
+   Int_t nt=0;
+   Double_t px,py,pz,x,y,z;
+   Int_t code,lab;
+   while (in>>lab>>code>>px>>py>>pz>>x>>y>>z) {
+      if (good[lab] != 0x3F) continue;
+      TParticle *p = (TParticle*)gAlice->Particle(lab);
+      gt[nt].lab=lab;
+      gt[nt].code=p->GetPdgCode();
+//**** px py pz - in global coordinate system
+      gt[nt].px=p->Px(); gt[nt].py=p->Py(); gt[nt].pz=p->Pz();
+      gt[nt].x=gt[nt].y=gt[nt].z=0.;
+      nt++;
+      if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
+   }
+
+   delete[] good;
+
+   delete gAlice; gAlice=0;
+   file->Close();
+
+   return nt;
+}
+
+
diff --git a/ITS/AliITSFindClustersV2.C b/ITS/AliITSFindClustersV2.C
new file mode 100644 (file)
index 0000000..46395e1
--- /dev/null
@@ -0,0 +1,172 @@
+#ifndef __CINT__
+  #include "AliRun.h"
+  #include "AliITS.h"
+  #include "AliITSgeom.h"
+  #include "AliITSRecPoint.h"
+  #include "AliITSclusterV2.h"
+  #include "AliITSsimulationFastPoints.h"
+
+  #include "TFile.h"
+  #include "TTree.h"
+  #include "TParticle.h"
+#endif
+
+Int_t AliITSFindClustersV2() {
+/****************************************************************
+ *  Just something to start with                                *
+ ****************************************************************/
+   cerr<<"Looking for clusters...\n";
+
+   if (gAlice) {delete gAlice; gAlice=0;}
+
+   TFile *in=TFile::Open("galice.root","update");
+   if (!in->IsOpen()) {
+      cerr<<"Can't open galice.root !\n"; 
+      return 1;
+   }
+   
+   if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
+      cerr<<"Can't find gAlice !\n";
+      return 2;
+   }
+
+   Int_t ev=0;
+   gAlice->GetEvent(ev);
+
+   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
+   if (!ITS) {
+      cerr<<"Can't find the ITS !\n";
+      return 3;
+   }
+
+   gAlice->MakeTree("R");
+   ITS->MakeBranch("R",0);
+   gAlice->TreeR()->Fill();
+
+   //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
+   AliITSsimulationFastPoints *sim = new AliITSsimulationFastPoints();
+   for (Int_t i=0;i<3;i++) { ITS->SetSimulationModel(i,sim); }
+
+   Int_t nsignal=25;
+   Int_t size=-1;
+   Int_t bgr_ev=Int_t(ev/nsignal);
+   ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
+   //////////////////////////////////////////////////////////////////////////
+
+   delete gAlice; gAlice=0;
+   in->Close();
+
+///////////////// Conversion AliITSRecPoint -> AliITSclusterV2 //////////////
+   /*TFile */in=TFile::Open("galice.root");
+
+   if (gAlice) {delete gAlice; gAlice=0;}
+
+   if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
+      cerr<<"Can't find gAlice !\n";
+      return 4;
+   }
+
+   gAlice->GetEvent(0);
+
+   /*AliITS */ITS  = (AliITS*)gAlice->GetModule("ITS");
+   if (!ITS) {
+      cerr<<"Can't find the ITS !\n";
+      return 5;
+   }
+   AliITSgeom *geom=ITS->GetITSgeom();
+   TFile *out=TFile::Open("AliITSclustersV2.root","new");
+   if (!out->IsOpen()) {
+      cerr<<"Delete old AliITSclustersV2.root !\n"; 
+      return 6;
+   }
+   geom->Write();
+
+   TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
+   TTree *cTree=new TTree("cTree","ITS clusters");
+   cTree->Branch("Clusters",&clusters);
+
+   TTree *pTree=gAlice->TreeR();
+   if (!pTree) { 
+      cerr<<"Can't get TreeR !\n";
+      return 7;
+   }
+   TBranch *branch=pTree->GetBranch("ITSRecPoints");
+   if (!branch) { 
+      cerr<<"Can't get ITSRecPoints branch !\n";
+      return 8;
+   }
+   TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
+   branch->SetAddress(&points);
+
+   TClonesArray &cl=*clusters;
+   Int_t nclusters=0;
+   Int_t nentr=(Int_t)pTree->GetEntries();
+
+   cerr<<"Number of entries: "<<nentr<<endl;
+
+   for (Int_t i=0; i<nentr; i++) {
+       if (!pTree->GetEvent(i)) {cTree->Fill(); continue;}
+       Int_t lay,lad,det; geom->GetModuleId(i-1,lay,lad,det);
+       Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); 
+       Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
+       Double_t yshift = x*rot[0] + y*rot[1];
+       Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
+       Int_t ncl=points->GetEntriesFast();
+       nclusters+=ncl;
+       for (Int_t j=0; j<ncl; j++) {
+          AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
+          Float_t lp[5];
+          lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
+          lp[1]=p->GetZ()+zshift;
+          lp[2]=p->GetSigmaX2();
+          lp[3]=p->GetSigmaZ2();
+          lp[4]=p->GetQ();
+          Int_t lab[6]; 
+          lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
+          lab[3]=ndet;
+
+          Int_t label=lab[0];
+          TParticle *part=(TParticle*)gAlice->Particle(label);
+          label=-3;
+          while (part->P() < 0.005) {
+             Int_t m=part->GetFirstMother();
+             if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
+             label=m;
+             part=(TParticle*)gAlice->Particle(label);
+          }
+          if      (lab[1]<0) lab[1]=label;
+          else if (lab[2]<0) lab[2]=label;
+          else cerr<<"No empty labels !\n";
+
+          new(cl[j]) AliITSclusterV2(lab,lp);
+       }
+       cTree->Fill(); clusters->Delete();
+       points->Delete();
+   }
+   cTree->Write();
+
+   cerr<<"Number of clusters: "<<nclusters<<endl;
+
+   delete cTree; delete clusters; delete points;
+
+   delete gAlice; gAlice=0;
+   in->Close();
+   out->Close();
+
+   return 0;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ITS/AliITSFindTracksV2.C b/ITS/AliITSFindTracksV2.C
new file mode 100644 (file)
index 0000000..c18e538
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef __CINT__
+  #include "AliITStrackerV2.h"
+
+  #include "TFile.h"
+  #include "TStopwatch.h"
+#endif
+
+Int_t AliITSFindTracksV2() {
+   cerr<<"Looking for tracks...\n";
+
+   TFile *out=TFile::Open("AliITStracksV2.root","new");
+   if (!out->IsOpen()) {cerr<<"Delete old AliITStracksV2.root !\n"; return 1;}
+
+   TFile *in=TFile::Open("AliTPCtracks.root");
+   if (!in->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 2;}
+
+   TFile *file=TFile::Open("AliITSclustersV2.root");
+   if (!file->IsOpen()) {cerr<<"Can't open AliITSclustersV2.root !\n";return 3;}
+
+   AliITSgeom *geom=(AliITSgeom*)file->Get("AliITSgeom");
+
+   TStopwatch timer;
+   AliITStrackerV2 tracker(geom);
+   Int_t rc=tracker.Clusters2Tracks(in,out);
+   timer.Stop(); timer.Print();
+
+   file->Close();
+   in->Close();
+   out->Close();
+
+   return rc;
+}
diff --git a/ITS/AliITSclusterV2.cxx b/ITS/AliITSclusterV2.cxx
new file mode 100644 (file)
index 0000000..1d1b43a
--- /dev/null
@@ -0,0 +1,24 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+//              Implementation of the ITS cluster class
+//
+//         Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+//-------------------------------------------------------------------------
+
+#include "AliITSclusterV2.h"
+
+ClassImp(AliITSclusterV2)
diff --git a/ITS/AliITSclusterV2.h b/ITS/AliITSclusterV2.h
new file mode 100644 (file)
index 0000000..e49e4d6
--- /dev/null
@@ -0,0 +1,39 @@
+#ifndef ALIITSCLUSTERV2_H
+#define ALIITSCLUSTERV2_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                        ITS Cluster Class
+//
+//      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
+//-------------------------------------------------------------------------
+
+#include "AliCluster.h"
+#include "AliITSrecoV2.h"
+
+//_____________________________________________________________________________
+class AliITSclusterV2 : public AliCluster {
+public:
+  AliITSclusterV2() : AliCluster() {fQ=0;}
+  AliITSclusterV2(Int_t *lab,Float_t *hit) : AliCluster(lab,hit) {
+    fIndex=lab[3];
+    fQ=hit[4];
+  }
+  void Use() {fSigmaY2=-fSigmaY2;}
+  void SetQ(Float_t q) {fQ=q;}
+
+  Int_t IsUsed() const {return (fSigmaY2<0) ? 1 : 0;}
+  Float_t GetQ() const {return fQ;}
+  Int_t GetDetectorIndex() const { return fIndex; }
+
+private:
+  Int_t    fIndex;    // detector index
+  Float_t  fQ ;       // Q of cluster (in ADC counts)
+  
+  ClassDef(AliITSclusterV2,1)  // ITS clusters
+};
+
+#endif
+
+
diff --git a/ITS/AliITSrecoV2.h b/ITS/AliITSrecoV2.h
new file mode 100644 (file)
index 0000000..a171b9c
--- /dev/null
@@ -0,0 +1,40 @@
+#ifndef ALIITSRECO_H
+#define ALIITSRECO_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                   ITS reconstruction name space
+//
+//       Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
+//-------------------------------------------------------------------------
+#include <Rtypes.h>
+#include <iostream.h>
+
+//namespace AliITSreco {    
+   const Int_t kMaxClusterPerLayer=3500*10;
+   const Int_t kMaxDetectorPerLayer=1000;
+   const Int_t kLayersToSkip=0;
+
+   const Int_t kMaxLayer=6;
+   const Double_t kSigmaY2[kMaxLayer]={
+      1.44e-6, 1.44e-6, 1.444e-5, 1.444e-5, 4.0e-6, 4.0e-6 
+   };
+   const Double_t kSigmaZ2[kMaxLayer]={
+     //4.9e-5, 4.9e-5, 7.84e-6, 7.84e-6, 0.006889, 0.006889
+     1.44e-4, 1.44e-4, 7.84e-6, 7.84e-6, 0.006889, 0.006889
+   };
+
+   const Double_t kChi2PerCluster=7.;//10.;//7
+   const Double_t kMaxChi2=15.;//20.; //15.
+   const Double_t kMaxRoad=13.;
+
+   const Double_t kSigmaYV=0.005e+0;
+   const Double_t kSigmaZV=0.010e+0;
+
+   const Double_t kConvConst=100/0.299792458/0.2; 
+//}
+
+//using namespace AliITSreco;
+
+#endif
diff --git a/ITS/AliITStestV2.C b/ITS/AliITStestV2.C
new file mode 100644 (file)
index 0000000..197a120
--- /dev/null
@@ -0,0 +1,14 @@
+Int_t AliITStestV2() {
+   Int_t rc=0;
+
+   gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersV2.C");
+   if (rc=AliITSFindClustersV2()) return rc;
+
+   gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindTracksV2.C");
+   if (rc=AliITSFindTracksV2()) return rc;
+
+   gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSComparisonV2.C");
+   if (rc=AliITSComparisonV2()) return rc;
+
+   return rc;
+}
diff --git a/ITS/AliITStrackV2.cxx b/ITS/AliITStrackV2.cxx
new file mode 100644 (file)
index 0000000..4dda61c
--- /dev/null
@@ -0,0 +1,800 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+//                Implementation of the ITS track class
+//
+//          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+//-------------------------------------------------------------------------
+
+#include <TMatrixD.h>
+
+#include <TMath.h>
+#include <iostream.h>
+
+#include "AliCluster.h"
+#include "../TPC/AliTPCtrack.h"
+#include "AliITStrackV2.h"
+
+ClassImp(AliITStrackV2)
+
+const Int_t kWARN=1;
+
+//____________________________________________________________________________
+AliITStrackV2::AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *) {
+  //------------------------------------------------------------------
+  //Convertion TPC track -> ITS track
+  //------------------------------------------------------------------
+  SetLabel(t.GetLabel());
+  SetChi2(0.);
+  SetNumberOfClusters(0);
+  fdEdx  = 0.;
+  fAlpha = t.GetAlpha();
+  if      (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
+  else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
+
+  //Convertion of the track parameters
+  Double_t x,p[5]; t.GetExternalParameters(x,p);
+  fX=x;    x=kConversionConstant;
+  fP0=p[0]; 
+  fP1=p[1]; 
+  fP2=p[2];
+  fP3=p[3];
+  fP4=p[4]/x; 
+
+  //Convertion of the covariance matrix
+  Double_t c[15]; t.GetExternalCovariance(c);
+
+  fC00=c[0 ];
+  fC10=c[1 ];   fC11=c[2 ];
+  fC20=c[3 ];   fC21=c[4 ];   fC22=c[5 ];
+  fC30=c[6 ];   fC31=c[7 ];   fC32=c[8 ];   fC33=c[9 ];
+  fC40=c[10]/x; fC41=c[11]/x; fC42=c[12]/x; fC43=c[13]/x; fC44=c[14]/x/x;
+
+  if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
+}
+
+//____________________________________________________________________________
+AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : AliKalmanTrack(t) {
+  //------------------------------------------------------------------
+  //Copy constructor
+  //------------------------------------------------------------------
+  fX=t.fX;
+  fAlpha=t.fAlpha;
+  fdEdx=t.fdEdx;
+
+  fP0=t.fP0; fP1=t.fP1; fP2=t.fP2; fP3=t.fP3; fP4=t.fP4;
+
+  fC00=t.fC00;
+  fC10=t.fC10;  fC11=t.fC11;
+  fC20=t.fC20;  fC21=t.fC21;  fC22=t.fC22;
+  fC30=t.fC30;  fC31=t.fC31;  fC32=t.fC32;  fC33=t.fC33;
+  fC40=t.fC40;  fC41=t.fC41;  fC42=t.fC42;  fC43=t.fC43;  fC44=t.fC44;
+
+  Int_t n=GetNumberOfClusters();
+  for (Int_t i=0; i<n; i++) fIndex[i]=t.fIndex[i];
+}
+
+//_____________________________________________________________________________
+Int_t AliITStrackV2::Compare(const TObject *o) const {
+  //-----------------------------------------------------------------
+  // This function compares tracks according to the their curvature
+  //-----------------------------------------------------------------
+  AliTPCtrack *t=(AliTPCtrack*)o;
+  Double_t co=TMath::Abs(t->Get1Pt());
+  Double_t c =TMath::Abs(Get1Pt());
+  if (c>co) return 1;
+  else if (c<co) return -1;
+  return 0;
+}
+
+//_____________________________________________________________________________
+void AliITStrackV2::GetExternalCovariance(Double_t cc[15]) const {
+  //-------------------------------------------------------------------------
+  // This function returns an external representation of the covriance matrix.
+  //   (See comments in AliTPCtrack.h about external track representation)
+  //-------------------------------------------------------------------------
+  Double_t a=kConvConst;
+
+  cc[0 ]=fC00;
+  cc[1 ]=fC10;   cc[2 ]=fC11;
+  cc[3 ]=fC20;   cc[4 ]=fC21;   cc[5 ]=fC22;
+  cc[6 ]=fC30;   cc[7 ]=fC31;   cc[8 ]=fC32;   cc[9 ]=fC33;
+  cc[10]=fC40*a; cc[11]=fC41*a; cc[12]=fC42*a; cc[13]=fC43*a; cc[14]=fC44*a*a;
+}
+
+//____________________________________________________________________________
+Int_t AliITStrackV2::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm) {
+  //------------------------------------------------------------------
+  //This function propagates a track to the minimal distance from the origin
+  //------------------------------------------------------------------
+  Double_t xv=fP2*(fX*fP2 - fP0*TMath::Sqrt(1.- fP2*fP2)); //linear approxim.
+  Propagate(fAlpha,xv,0.,0.,pm);   
+  return 0;
+}
+
+//____________________________________________________________________________
+Int_t AliITStrackV2::
+GetGlobalXYZat(Double_t xk, Double_t &x, Double_t &y, Double_t &z) const {
+  //------------------------------------------------------------------
+  //This function returns a track position in the global system
+  //------------------------------------------------------------------
+  Double_t dx=xk-fX;
+  Double_t f1=fP2, f2=f1 + fP4*dx;
+  if (TMath::Abs(f2) >= 0.99999) {
+     Int_t n=GetNumberOfClusters();
+    if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Propagation failed !\n";
+    return 0;
+  }
+
+  Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
+  
+  Double_t yk = fP0 + dx*(f1+f2)/(r1+r2);
+  Double_t zk = fP1 + dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
+
+  Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
+  x = xk*cs - yk*sn;
+  y = xk*sn + yk*cs;
+  z = zk;
+
+  return 1;
+}
+
+//_____________________________________________________________________________
+Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const 
+{
+  //-----------------------------------------------------------------
+  // This function calculates a predicted chi2 increment.
+  //-----------------------------------------------------------------
+  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
+  r00+=fC00; r01+=fC10; r11+=fC11;
+
+  Double_t det=r00*r11 - r01*r01;
+  if (TMath::Abs(det) < 1.e-10) {
+    Int_t n=GetNumberOfClusters();
+    if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
+    return 1e10;
+  }
+  Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
+  
+  Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
+  
+  return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
+}
+
+//_____________________________________________________________________________
+Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c,Double_t *m,
+Double_t x0, Double_t pm=0.139) const {
+  //-----------------------------------------------------------------
+  // This function calculates a chi2 increment with a vertex contraint 
+  //-----------------------------------------------------------------
+  TVectorD x(5); x(0)=fP0; x(1)=fP1; x(2)=fP2; x(3)=fP3; x(4)=fP4;
+  TMatrixD C(5,5);
+  C(0,0)=fC00; 
+  C(1,0)=fC10; C(1,1)=fC11; 
+  C(2,0)=fC20; C(2,1)=fC21; C(2,2)=fC22;
+  C(3,0)=fC30; C(3,1)=fC31; C(3,2)=fC32; C(3,3)=fC33;
+  C(4,0)=fC40; C(4,1)=fC41; C(4,2)=fC42; C(4,3)=fC43; C(4,4)=fC44;
+
+  C(0,1)=C(1,0);
+  C(0,2)=C(2,0); C(1,2)=C(2,1);
+  C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
+  C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
+
+  TMatrixD H(4,5); H.UnitMatrix();
+  Double_t dy=(c->GetY() - m[0]), dz=(c->GetZ() - m[1]);
+
+  Double_t dr=TMath::Sqrt(fX*fX + dy*dy);
+  Double_t r =TMath::Sqrt(4/dr/dr - fP4*fP4);
+  Double_t sn=0.5*(fP4*fX + dy*r);
+  Double_t tg=0.5*fP4*dz/TMath::ASin(0.5*fP4*dr);
+  TVectorD mm(4); 
+  mm(0)=m[0]=c->GetY(); mm(1)=m[1]=c->GetZ(); mm(2)=m[2]=sn; mm(3)=m[3]=tg;
+
+  Double_t v22=0.,v33=0.;
+  //x0=0.;
+  if (x0!=0.) {
+     Double_t pp2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
+     Double_t beta2=pp2/(pp2 + pm*pm);
+     x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
+     Double_t theta2=14.1*14.1/(beta2*pp2*1e6)*x0;
+     v22 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
+     v33 = theta2*(1.+ GetTgl()*GetTgl())*(1. + GetTgl()*GetTgl());
+  }
+  Double_t sy2=c->GetSigmaY2(), sz2=c->GetSigmaZ2();
+  v22+=kSigmaYV*kSigmaYV/dr/dr;
+  v22+=sy2/dr/dr;
+  Double_t v20=sy2/dr;
+
+  v33+=kSigmaZV*kSigmaZV/dr/dr;
+  v33+=sz2/dr/dr;
+  Double_t v31=sz2/dr;
+
+  TMatrixD V(4,4); 
+  V(0,0)=m[4 ]=sy2; V(0,1)=m[5 ]=0.;  V(0,2)=m[6 ]=v20; V(0,3)=m[7 ]=0.;
+  V(1,0)=m[8 ]=0.;  V(1,1)=m[9 ]=sz2; V(1,2)=m[10]=0.;  V(1,3)=m[11]=v31;
+  V(2,0)=m[12]=v20; V(2,1)=m[13]=0.;  V(2,2)=m[14]=v22; V(2,3)=m[15]=0.;
+  V(3,0)=m[16]=0.;  V(3,1)=m[17]=v31; V(3,2)=m[18]=0.;  V(3,3)=m[19]=v33;
+
+  TVectorD res=x;  res*=H; res-=mm; //res*=-1; 
+  TMatrixD tmp(H,TMatrixD::kMult,C);
+  TMatrixD R(tmp,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed,H)); R+=V;
+  
+  Double_t det=R.Determinant();
+  if (TMath::Abs(det) < 1.e-25) {
+    Int_t n=GetNumberOfClusters();
+    if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Singular matrix !\n";
+    return 1e10;
+  }
+
+  R.Invert();
+
+  TVectorD rs=res;
+  res*=R;
+  return rs*res;
+}
+
+//____________________________________________________________________________
+Int_t 
+AliITStrackV2::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) {
+  //------------------------------------------------------------------
+  //This function propagates a track
+  //------------------------------------------------------------------
+  Double_t x1=fX, x2=xk, dx=x2-x1;
+  Double_t f1=fP2, f2=f1 + fP4*dx;
+  if (TMath::Abs(f2) >= 0.99999) {
+    Int_t n=GetNumberOfClusters();
+    if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Propagation failed !\n";
+    return 0;
+  }
+
+  Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
+  
+  fP0 += dx*(f1+f2)/(r1+r2);
+  fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
+  fP2 += dx*fP4;
+
+  //f = F - 1
+  
+  Double_t f02=    dx/(r1*r1*r1);
+  Double_t f04=0.5*dx*dx/(r1*r1*r1);
+  Double_t f12=    dx*fP3*f1/(r1*r1*r1);
+  Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1);
+  Double_t f13=    dx/r1;
+  Double_t f24=    dx; 
+  
+  //b = C*ft
+  Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
+  Double_t b02=f24*fC40;
+  Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
+  Double_t b12=f24*fC41;
+  Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
+  Double_t b22=f24*fC42;
+  Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
+  Double_t b42=f24*fC44;
+  Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
+  Double_t b32=f24*fC43;
+  
+  //a = f*b = f*C*ft
+  Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
+  Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
+  Double_t a22=f24*b42;
+
+  //F*C*Ft = C + (b + bt + a)
+  fC00 += b00 + b00 + a00;
+  fC10 += b10 + b01 + a01; 
+  fC20 += b20 + b02 + a02;
+  fC30 += b30;
+  fC40 += b40;
+  fC11 += b11 + b11 + a11;
+  fC21 += b21 + b12 + a12;
+  fC31 += b31; 
+  fC41 += b41;
+  fC22 += b22 + b22 + a22;
+  fC32 += b32;
+  fC42 += b42;
+
+  fX=x2;
+
+  Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
+  Double_t beta2=p2/(p2 + pm*pm);
+
+  //Multiple scattering******************
+  //x0=0.;
+  if (x0!=0) {
+     x0*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
+     Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
+     fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
+     fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
+     fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
+     fC44 += theta2*fP3*fP4*fP3*fP4;
+  }
+
+  //Energy losses************************
+  if (rho!=0.) {
+     rho*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
+     Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*rho;
+     if (x1 < x2) dE=-dE;
+     fP4*=(1.- sqrt(p2+pm*pm)/p2*dE);
+  }
+
+  if (!Invariant()) {cout<<"Propagate !\n"; return 0;}
+
+  return 1;
+}
+
+//____________________________________________________________________________
+Int_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, UInt_t index) {
+  //------------------------------------------------------------------
+  //This function updates track parameters
+  //------------------------------------------------------------------
+  Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
+  Double_t c00=fC00;
+  Double_t c10=fC10, c11=fC11;
+  Double_t c20=fC20, c21=fC21, c22=fC22;
+  Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
+  Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
+
+
+  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
+  r00+=fC00; r01+=fC10; r11+=fC11;
+  Double_t det=r00*r11 - r01*r01;
+  Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
+
+  Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
+  Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
+  Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
+  Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
+  Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
+
+  Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
+  Double_t sf=fP2 + k20*dy + k21*dz;
+  /*
+  if (TMath::Abs(sf) >= 0.99999) {
+    Int_t n=GetNumberOfClusters();
+    if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Filtering failed !\n";
+    return 0;
+  }
+  */
+  fP0 += k00*dy + k01*dz;
+  fP1 += k10*dy + k11*dz;
+  fP2  = sf;
+  fP3 += k30*dy + k31*dz;
+  fP4 += k40*dy + k41*dz;
+  
+  Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
+  Double_t c12=fC21, c13=fC31, c14=fC41;
+
+  fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
+  fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
+  fC40-=k00*c04+k01*c14; 
+
+  fC11-=k10*c01+k11*fC11;
+  fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
+  fC41-=k10*c04+k11*c14; 
+
+  fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
+  fC42-=k20*c04+k21*c14; 
+
+  fC33-=k30*c03+k31*c13;
+  fC43-=k30*c04+k31*c14; 
+
+  fC44-=k40*c04+k41*c14; 
+
+  if (!Invariant()) {
+     fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
+     fC00=c00;
+     fC10=c10; fC11=c11;
+     fC20=c20; fC21=c21; fC22=c22;
+     fC30=c30; fC31=c31; fC32=c32; fC33=c33;
+     fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
+     return 0;
+  }
+
+  Int_t n=GetNumberOfClusters();
+  fIndex[n]=index;
+  SetNumberOfClusters(n+1);
+  SetChi2(GetChi2()+chi2);
+
+  return 1;
+}
+
+
+//____________________________________________________________________________
+Int_t AliITStrackV2::Update(const Double_t* m, Double_t chi2, UInt_t index) {
+  //------------------------------------------------------------------
+  //This function updates track parameters with a vertex constraint
+  //------------------------------------------------------------------
+  Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
+  Double_t c00=fC00;
+  Double_t c10=fC10, c11=fC11;
+  Double_t c20=fC20, c21=fC21, c22=fC22;
+  Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
+  Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
+
+
+  TVectorD x(5); x(0)=fP0; x(1)=fP1; x(2)=fP2; x(3)=fP3; x(4)=fP4;
+  TMatrixD C(5,5);
+  C(0,0)=fC00; 
+  C(1,0)=fC10; C(1,1)=fC11; 
+  C(2,0)=fC20; C(2,1)=fC21; C(2,2)=fC22;
+  C(3,0)=fC30; C(3,1)=fC31; C(3,2)=fC32; C(3,3)=fC33;
+  C(4,0)=fC40; C(4,1)=fC41; C(4,2)=fC42; C(4,3)=fC43; C(4,4)=fC44;
+
+  C(0,1)=C(1,0);
+  C(0,2)=C(2,0); C(1,2)=C(2,1);
+  C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
+  C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
+
+  TMatrixD H(4,5); H.UnitMatrix();
+  TMatrixD Ht(TMatrixD::kTransposed,H);
+  TVectorD mm(4); mm(0)=m[0]; mm(1)=m[1]; mm(2)=m[2]; mm(3)=m[3];
+  TMatrixD V(4,4); 
+  V(0,0)=m[4 ]; V(0,1)=m[5 ]; V(0,2)=m[6 ]; V(0,3)=m[7 ];
+  V(1,0)=m[8 ]; V(1,1)=m[9 ]; V(1,2)=m[10]; V(1,3)=m[11];
+  V(2,0)=m[12]; V(2,1)=m[13]; V(2,2)=m[14]; V(2,3)=m[15];
+  V(3,0)=m[16]; V(3,1)=m[17]; V(3,2)=m[18]; V(3,3)=m[19];
+
+  TMatrixD tmp(H,TMatrixD::kMult,C);
+  TMatrixD R(tmp,TMatrixD::kMult,Ht); R+=V;
+
+  R.Invert();
+  
+  TMatrixD K(C,TMatrixD::kMult,Ht); K*=R;
+  
+  TVectorD savex=x;
+  x*=H; x-=mm; x*=-1; x*=K; x+=savex;
+
+  TMatrixD saveC=C;
+  C.Mult(K,tmp); C-=saveC; C*=-1;
+
+  fP0=x(0); fP1=x(1); fP2=x(2); fP3=x(3); fP4=x(4);
+  fC00=C(0,0); 
+  fC10=C(1,0); fC11=C(1,1); 
+  fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
+  fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
+  fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
+
+
+  if (!Invariant()) {
+     fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
+     fC00=c00;
+     fC10=c10; fC11=c11;
+     fC20=c20; fC21=c21; fC22=c22;
+     fC30=c30; fC31=c31; fC32=c32; fC33=c33;
+     fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
+     return 0;
+  }
+
+  Int_t n=GetNumberOfClusters();
+  fIndex[n]=index;
+  SetNumberOfClusters(n+1);
+  SetChi2(GetChi2()+chi2);
+
+  return 1;
+}
+
+Int_t AliITStrackV2::Invariant() const {
+  //------------------------------------------------------------------
+  // This function is for debugging purpose only
+  //------------------------------------------------------------------
+  //if (TMath::Abs(fP1)>11.5)
+  //if (fP1*fP4<0) {cout<<"fP1*fP4="<<fP1*fP4<<' '<<fP1<<endl; return 0;}
+  if (TMath::Abs(fP2)>=1) {cout<<"fP2="<<fP2<<endl; return 0;}
+
+  if (fC00<=0) {cout<<"fC00="<<fC00<<endl; return 0;}
+  if (fC11<=0) {cout<<"fC11="<<fC11<<endl; return 0;}
+  if (fC22<=0) {cout<<"fC22="<<fC22<<endl; return 0;}
+  if (fC33<=0) {cout<<"fC33="<<fC33<<endl; return 0;}
+  if (fC44<=0) {cout<<"fC44="<<fC44<<endl; return 0;}
+
+   TMatrixD m(5,5);
+   m(0,0)=fC00; 
+   m(1,0)=fC10; m(1,1)=fC11; 
+   m(2,0)=fC20; m(2,1)=fC21; m(2,2)=fC22;
+   m(3,0)=fC30; m(3,1)=fC31; m(3,2)=fC32; m(3,3)=fC33;
+   m(4,0)=fC40; m(4,1)=fC41; m(4,2)=fC42; m(4,3)=fC43; m(4,4)=fC44;
+
+   m(0,1)=m(1,0);
+   m(0,2)=m(2,0); m(1,2)=m(2,1);
+   m(0,3)=m(3,0); m(1,3)=m(3,1); m(2,3)=m(3,2);
+   m(0,4)=m(4,0); m(1,4)=m(4,1); m(2,4)=m(4,2); m(3,4)=m(4,3);
+   /*
+   Double_t det=m.Determinant(); 
+
+   if (det <= 0) {
+       cout<<" bad determinant "<<det<<endl;
+       m.Print(); 
+       return 0;
+   }
+   */
+   return 1;
+}
+
+//____________________________________________________________________________
+Int_t AliITStrackV2::Propagate(Double_t alp, Double_t xk,
+Double_t x0,Double_t rho,Double_t pm) {
+  //------------------------------------------------------------------
+  //This function propagates a track
+  //------------------------------------------------------------------
+  Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
+  Double_t c00=fC00;
+  Double_t c10=fC10, c11=fC11;
+  Double_t c20=fC20, c21=fC21, c22=fC22;
+  Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
+  Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
+
+
+  Double_t dalp=alp-fAlpha;
+
+  Double_t ca=TMath::Cos(dalp), sa=TMath::Sin(dalp);
+  Double_t sf=fP2, cf=TMath::Sqrt(1.- fP2*fP2);  
+
+  Double_t pp2=fP2*ca - cf*sa;
+  if (TMath::Abs(pp2) >= 0.99999) {
+     Int_t n=GetNumberOfClusters();
+     if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Rotation failed !\n";
+     return 0;
+  }
+
+  fAlpha = alp;
+  if      (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
+  else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
+  
+  Double_t x1=fX, y1=fP0;
+
+  fX = x1*ca + y1*sa;
+  fP0=-x1*sa + y1*ca;
+  fP2 = pp2;
+
+  cf=ca + sf*sa/cf;
+
+  if (!Invariant()) {cout<<dalp<<" Rotate !\n"; return 0;}
+
+  x1=fX; Double_t x2=xk, dx=x2-x1;
+  Double_t f1=fP2, f2=f1 + fP4*dx;
+  if (TMath::Abs(f2) >= 0.99999) {
+    Int_t n=GetNumberOfClusters();
+    if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Propagation failed !\n";
+    return 0;
+  }
+
+  Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
+  
+  fP0 += dx*(f1+f2)/(r1+r2);
+  fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
+  fP2 += dx*fP4;
+
+  //f = F - 1
+  Double_t f02=    dx/(r1*r1*r1);
+  Double_t f04=0.5*dx*dx/(r1*r1*r1);
+  Double_t f12=    dx*fP3*f1/(r1*r1*r1);
+  Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1);
+  Double_t f13=    dx/r1;
+  Double_t f24=    dx; 
+  /*
+  //b = C*ft
+  Double_t b00=f02*fC20 + f03*fC30, b01=f12*fC20 + f13*fC30 + f14*fC40;
+  Double_t b02=f23*fC30;
+  Double_t b10=f02*fC21 + f03*fC31, b11=f12*fC21 + f13*fC31 + f14*fC41;
+  Double_t b12=f23*fC31;
+  Double_t b20=f02*fC22 + f03*fC32, b21=f12*fC22 + f13*fC32 + f14*fC42;
+  Double_t b22=f23*fC32;
+  Double_t b30=f02*fC32 + f03*fC33, b31=f12*fC32 + f13*fC33 + f14*fC43;
+  Double_t b32=f23*fC33;
+  Double_t b40=f02*fC42 + f03*fC43, b41=f12*fC42 + f13*fC43 + f14*fC44;
+  Double_t b42=f23*fC43;
+  
+  //a = f*b = f*C*ft
+  Double_t a00=f02*b20+f03*b30,a01=f02*b21+f03*b31,a02=f02*b22+f03*b32;
+  Double_t a11=f12*b21+f13*b31+f14*b41,a12=f12*b22+f13*b32+f14*b42;
+  Double_t a22=f23*b32;
+
+  //F*C*Ft = C + (b + bt + a)
+  fC00 += b00 + b00 + a00;
+  fC10 += b10 + b01 + a01; 
+  fC20 += b20 + b02 + a02;
+  fC30 += b30;
+  fC40 += b40;
+  fC11 += b11 + b11 + a11;
+  fC21 += b21 + b12 + a12;
+  fC31 += b31; 
+  fC41 += b41;
+  fC22 += b22 + b22 + a22;
+  fC32 += b32;
+  fC42 += b42;
+*/
+
+ TMatrixD F(5,5); F.UnitMatrix();
+ F(0,0)=-(f1+f2)/(r1+r2)*sa + ca; F(0,2)=f02*cf; F(0,4)=f04;
+ F(1,0)=-(f1+f2)/(f1*r2 + f2*r1)*fP3*sa; F(1,2)=f12*cf; F(1,4)=f14; F(1,3)=f13;
+ F(2,0)=-fP4*sa; F(2,2)=cf; F(2,4)=f24;
+
+  TMatrixD C(5,5);
+  C(0,0)=fC00; 
+  C(1,0)=fC10; C(1,1)=fC11; 
+  C(2,0)=fC20; C(2,1)=fC21; C(2,2)=fC22;
+  C(3,0)=fC30; C(3,1)=fC31; C(3,2)=fC32; C(3,3)=fC33;
+  C(4,0)=fC40; C(4,1)=fC41; C(4,2)=fC42; C(4,3)=fC43; C(4,4)=fC44;
+
+  C(0,1)=C(1,0);
+  C(0,2)=C(2,0); C(1,2)=C(2,1);
+  C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
+  C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
+
+  TMatrixD tmp(C,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
+  C.Mult(F,tmp);
+
+  fC00=C(0,0); 
+  fC10=C(1,0); fC11=C(1,1); 
+  fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
+  fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
+  fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
+
+  pp2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
+  Double_t beta2=pp2/(pp2 + pm*pm);
+
+  //Multiple scattering******************
+  //x0=0.;
+  if (x0!=0.) {
+     x0*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
+     Double_t theta2=14.1*14.1/(beta2*pp2*1e6)*x0;
+     fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
+     fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
+     fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
+     fC44 += theta2*fP3*fP4*fP3*fP4;
+  }
+
+  //Energy losses************************
+  if (rho!=0.) {  
+     rho*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
+     Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*rho;
+     if (x1 < x2) dE=-dE;
+     fP4*=(1.- sqrt(pp2+pm*pm)/pp2*dE);
+  }
+
+  if (!Invariant()) {
+     fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
+     fC00=c00;
+     fC10=c10; fC11=c11;
+     fC20=c20; fC21=c21; fC22=c22;
+     fC30=c30; fC31=c31; fC32=c32; fC33=c33;
+     fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
+     return 0;
+  }
+
+  fX=x2;
+
+  return 1;
+}
+
+Double_t AliITStrackV2::GetD() const {
+  //------------------------------------------------------------------
+  //This function calculates the transverse impact parameter
+  //------------------------------------------------------------------
+  Double_t sn=fP4*fX - fP2, cs=fP4*fP0 + TMath::Sqrt(1.- fP2*fP2);
+  Double_t a=2*(fX*fP2 - fP0*TMath::Sqrt(1.- fP2*fP2))-fP4*(fX*fX + fP0*fP0);
+  if (fP4<0) a=-a;
+  return a/(1 + TMath::Sqrt(sn*sn + cs*cs));
+} 
+
+
+Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
+  //------------------------------------------------------------------
+  //This function improves angular track parameters  
+  //------------------------------------------------------------------
+  Double_t dy=fP0-yv, dz=fP1-zv;
+  Double_t r2=fX*fX+dy*dy;
+  Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
+  Double_t beta2=p2/(p2 + 0.14*0.14);
+  x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
+  Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
+
+  Double_t par=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
+  Double_t sigma2 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
+  sigma2 += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
+  sigma2 += kSigmaYV*kSigmaYV/r2;
+  sigma2 += 0.25*fC44*fX*fX;
+  Double_t eps2=sigma2/(fC22+sigma2), eps=TMath::Sqrt(eps2);
+  if (10*r2*fC44<fC22) {
+     fP2 = eps2*fP2 + (1-eps2)*par;
+     fC22*=eps2; fC21*=eps; fC20*=eps; fC32*=eps; fC42*=eps;
+  }
+
+  par=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
+  sigma2=theta2;
+  sigma2 += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
+  sigma2 += kSigmaZV*kSigmaZV/r2;
+  eps2=sigma2/(fC33+sigma2); eps=TMath::Sqrt(eps2);
+  Double_t tgl=fP3;
+  fP3 = eps2*fP3 + (1-eps2)*par;
+  fC33*=eps2; fC32*=eps; fC31*=eps; fC30*=eps; fC43*=eps;
+
+  eps=TMath::Sqrt((1+fP3*fP3)/(1+tgl*tgl));
+  fP4*=eps;
+  fC44*=eps*eps; fC43*=eps;fC42*=eps; fC41*=eps; fC40*=eps;
+
+  if (!Invariant()) return 0;
+  return 1;
+} 
+
+/*
+Int_t AliITStrackV2::Improve(Double_t x0,Double_t xv,Double_t yv) {
+  //------------------------------------------------------------------
+  //This function improves angular track parameters  
+  //------------------------------------------------------------------
+  TMatrixD I(5,5);
+  TVectorD v(5); v(0)=fP0; v(1)=fP1; v(2)=fP2; v(3)=fP3; v(4)=fP4;
+
+  Double_t r2=fX*fX+fP0*fP0;
+  Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
+  Double_t beta2=p2/(p2 + 0.14*0.14);
+  x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
+  Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
+
+  v(2)=0.5*(fP4*fX + fP0*TMath::Sqrt(4/r2-fP4*fP4));
+  Double_t sigma2 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
+  sigma2 += fC00/r2*(1.- fP0*fP0/r2)*(1.- fP0*fP0/r2);
+  sigma2 += kSigmaYV*kSigmaYV/r2;
+  I(2,2)=1/sigma2;
+
+  v(3)=0.5*fP4*fP1/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
+  sigma2=theta2;
+  sigma2 += fC11/r2+fC00*fP0*fP0*fP1*fP1/(r2*r2*r2);
+  sigma2 += kSigmaZV*kSigmaZV/r2;
+  I(3,3)=1/sigma2;
+
+  Double_t tgl=fP3;
+
+  TVectorD x(5); x(0)=fP0; x(1)=fP1; x(2)=fP2; x(3)=fP3; x(4)=fP4;
+  TMatrixD C(5,5);
+  C(0,0)=fC00; 
+  C(1,0)=fC10; C(1,1)=fC11; 
+  C(2,0)=fC20; C(2,1)=fC21; C(2,2)=fC22;
+  C(3,0)=fC30; C(3,1)=fC31; C(3,2)=fC32; C(3,3)=fC33;
+  C(4,0)=fC40; C(4,1)=fC41; C(4,2)=fC42; C(4,3)=fC43; C(4,4)=fC44;
+
+  C(0,1)=C(1,0);
+  C(0,2)=C(2,0); C(1,2)=C(2,1);
+  C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
+  C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
+
+  TMatrixD tmp(I,TMatrixD::kMult,C),U(5,5); U.UnitMatrix();
+  U+=tmp;
+  U.Invert();
+  TMatrixD W1(U);
+  TMatrixD W2(tmp,TMatrixD::kMult,W1);
+
+  v*=W2; x*=W1; x+=v;
+
+  C*=W1;
+
+
+  fP0=x(0); fP1=x(1); fP2=x(2); fP3=x(3); fP4=x(4);
+  fC00=C(0,0); 
+  fC10=C(1,0); fC11=C(1,1); 
+  fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
+  fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
+  fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
+
+  eps=TMath::Sqrt((1+fP3*fP3)/(1+tgl*tgl));
+  fP4*=eps;
+  fC44*=eps*eps; fC43*=eps;fC42*=eps; fC41*=eps; fC40*=eps;
+
+  if (!Invariant()) return 0;
+  return 1;
+} 
+*/
+
diff --git a/ITS/AliITStrackV2.h b/ITS/AliITStrackV2.h
new file mode 100644 (file)
index 0000000..2defb7a
--- /dev/null
@@ -0,0 +1,124 @@
+#ifndef ALIITSTRACKV2_H
+#define ALIITSTRACKV2_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                       ITS Track Class
+//
+//        Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
+//-------------------------------------------------------------------------
+
+
+/*****************************************************************************
+ *                          December 18, 2000                                *
+ *  Internal view of the ITS track parametrisation as well as the order of   *
+ *           track parameters are subject for possible changes !             *
+ *  Use GetExternalParameters() and GetExternalCovariance() to access ITS    *
+ *      track information regardless of its internal representation.         *
+ * This formation is now fixed in the following way:                         *
+ *      external param0:   local Y-coordinate of a track (cm)                *
+ *      external param1:   local Z-coordinate of a track (cm)                *
+ *      external param2:   local sine of the track momentum azimuthal angle  *
+ *      external param3:   tangent of the track momentum dip angle           *
+ *      external param4:   1/pt (1/(GeV/c))                                  *
+ *****************************************************************************/
+
+#include <AliKalmanTrack.h>
+
+#include "AliITSrecoV2.h"
+
+class AliTPCtrack;
+
+//_____________________________________________________________________________
+class AliITStrackV2  : public AliKalmanTrack {
+public:
+  AliITStrackV2():AliKalmanTrack(){}
+  AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *);
+  AliITStrackV2(const AliITStrackV2& t);
+  Int_t 
+  PropagateToVertex(Double_t x0=36.66,Double_t rho=1.2e-3,Double_t pm=0.139);
+  void SetdEdx(Float_t dedx) {fdEdx=dedx;}
+  void SetDetectorIndex(Int_t i) {SetLabel(i);}
+  
+  void *operator new(size_t s,void *p) { return p; }
+  void *operator new(size_t s) { return ::operator new(s); }
+  
+  Int_t GetDetectorIndex() const {return GetLabel();}
+  Double_t GetX()    const {return fX;}
+  Double_t GetAlpha()const {return fAlpha;}
+  Float_t  GetdEdx() const {return fdEdx;}
+
+  Double_t GetY()    const {return fP0;}
+  Double_t GetZ()    const {return fP1;}
+  Double_t GetSnp()  const {return fP2;}
+  Double_t GetTgl()  const {return fP3;}
+  Double_t Get1Pt()  const {return fP4*kConvConst;}
+
+
+  Double_t GetD() const;
+
+
+  Double_t GetSigmaY2() const {return fC00;}
+  Double_t GetSigmaZ2() const {return fC11;}
+
+  Int_t Compare(const TObject *o) const;
+
+  void GetExternalParameters(Double_t& xr, Double_t x[5]) const ;
+  void GetExternalCovariance(Double_t cov[15]) const ;
+
+  Int_t GetClusterIndex(Int_t i) const {return fIndex[i];}
+  Int_t GetGlobalXYZat(Double_t r,Double_t &x,Double_t &y,Double_t &z) const;
+
+  Int_t Propagate(Double_t alpha,
+                  Double_t xr,Double_t x0,Double_t rho,Double_t pm=0.139);
+
+  Double_t GetPredictedChi2(const AliCluster *cluster) const;
+  Int_t Update(const AliCluster* cl,Double_t chi2,UInt_t i);
+
+  Double_t GetPredictedChi2(const AliCluster *cluster, Double_t *m,
+                  Double_t x0, Double_t pm=0.139) const;
+  Int_t Update(const Double_t *m, Double_t chi2, UInt_t i);
+  Int_t Improve(Double_t x0,Double_t yv,Double_t zv);
+
+  Int_t Invariant() const;
+
+  //protected:
+Int_t 
+PropagateTo(Double_t xr,Double_t x0=21.82,Double_t rho=2.33,Double_t pm=0.139);
+private:
+  Double_t fX;              // X-coordinate of this track (reference plane)
+  Double_t fAlpha;          // rotation angle
+
+  Double_t fdEdx;           // dE/dx
+
+  Double_t fP0;             // Y-coordinate of a track 
+  Double_t fP1;             // Z-coordinate of a track
+  Double_t fP2;             // sine of the track momentum azimuthal angle
+  Double_t fP3;             // tangent of the track momentum dip angle
+  Double_t fP4;             // track curvature
+
+  Double_t fC00;                         // covariance
+  Double_t fC10, fC11;                   // matrix
+  Double_t fC20, fC21, fC22;             // of the
+  Double_t fC30, fC31, fC32, fC33;       // track
+  Double_t fC40, fC41, fC42, fC43, fC44; // parameters 
+
+  UInt_t fIndex[kMaxLayer]; // indices of associated clusters 
+
+  ClassDef(AliITStrackV2,1)   //ITS reconstructed track
+};
+
+inline 
+void AliITStrackV2::GetExternalParameters(Double_t& xr, Double_t x[5]) const {
+  //---------------------------------------------------------------------
+  // This function return external TPC track representation
+  //---------------------------------------------------------------------
+     xr=fX;          
+     x[0]=GetY(); x[1]=GetZ(); x[2]=GetSnp(); x[3]=GetTgl(); x[4]=Get1Pt();
+}
+
+#endif
+
+
diff --git a/ITS/AliITStrackerV2.cxx b/ITS/AliITStrackerV2.cxx
new file mode 100644 (file)
index 0000000..39b8e50
--- /dev/null
@@ -0,0 +1,620 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+//               Implementation of the ITS tracker class
+//
+//          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+//-------------------------------------------------------------------------
+#include <TFile.h>
+#include <TTree.h>
+#include <TRandom.h>
+#include <iostream.h>
+
+#include "AliITSgeom.h"
+#include "AliITSRecPoint.h"
+#include "../TPC/AliTPCtrack.h"
+#include "AliITSclusterV2.h"
+#include "AliITStrackerV2.h"
+
+//#define DEBUG
+
+#ifdef DEBUG
+Int_t LAB=236;
+#endif
+
+extern TRandom *gRandom;
+
+AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers
+
+AliITStrackerV2::
+AliITStrackerV2(const AliITSgeom *geom) throw (const Char_t *) {
+  //--------------------------------------------------------------------
+  //This is the AliITStracker constructor
+  //It also reads clusters from a root file
+  //--------------------------------------------------------------------
+  fYV=fZV=0.;
+
+  AliITSgeom *g=(AliITSgeom*)geom;
+
+  Float_t x,y,z;
+  Int_t i;
+  for (i=1; i<kMaxLayer+1; i++) {
+    Int_t nlad=g->GetNladders(i);
+    Int_t ndet=g->GetNdetectors(i);
+
+    g->GetTrans(i,1,1,x,y,z); 
+    Double_t r=TMath::Sqrt(x*x + y*y);
+    Double_t poff=TMath::ATan2(y,x);
+    Double_t zoff=z;
+
+    g->GetTrans(i,1,2,x,y,z);
+    r += TMath::Sqrt(x*x + y*y);
+    g->GetTrans(i,2,1,x,y,z);
+    r += TMath::Sqrt(x*x + y*y);
+    g->GetTrans(i,2,2,x,y,z);
+    r += TMath::Sqrt(x*x + y*y);
+    r*=0.25;
+
+    new (fLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
+
+    for (Int_t j=1; j<nlad+1; j++) {
+      for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
+        Float_t x,y,zshift; g->GetTrans(i,j,k,x,y,zshift); 
+        Double_t rot[9]; g->GetRotMatrix(i,j,k,rot);
+
+        Double_t r =-x*rot[1] + y*rot[0];         if (i==1) r=-r;
+        Double_t phi=TMath::ATan2(rot[1],rot[0]); if (i==1) phi-=3.1415927;
+        phi+=0.5*TMath::Pi(); 
+        AliITSdetector &det=fLayers[i-1].GetDetector((j-1)*ndet + k-1); 
+
+if (phi<0) phi += 2*TMath::Pi();
+
+        new(&det) AliITSdetector(r,phi); 
+      } 
+    }  
+
+  }
+
+  try {
+     //Read clusters
+     TTree *cTree=(TTree*)gDirectory->Get("cTree");
+     if (!cTree) throw 
+        ("AliITStrackerV2::AliITStrackerV2 can't get cTree !\n");
+
+     TBranch *branch=cTree->GetBranch("Clusters");
+     if (!branch) throw 
+        ("AliITStrackerV2::AliITStrackerV2 can't get Clusters branch !\n");
+
+     TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
+     branch->SetAddress(&clusters);
+
+     Int_t nentr=(Int_t)cTree->GetEntries();
+     for (i=0; i<nentr; i++) {
+       if (!cTree->GetEvent(i)) continue;
+       Int_t lay,lad,det; g->GetModuleId(i-1,lay,lad,det);
+       Int_t ncl=clusters->GetEntriesFast();
+       while (ncl--) {
+         AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
+
+#ifdef DEBUG
+if (c->GetLabel(2)!=LAB)
+if (c->GetLabel(1)!=LAB)
+if (c->GetLabel(0)!=LAB) continue;
+cout<<lay-1<<' '<<lad-1<<' '<<det-1<<' '<<c->GetY()<<' '<<c->GetZ()<<endl;
+#endif
+
+         fLayers[lay-1].InsertCluster(new AliITSclusterV2(*c));
+       }
+       clusters->Delete();
+     }
+  }
+  catch (const Char_t *msg) {
+    cerr<<msg<<endl;
+    throw;
+  }
+
+  fI=kMaxLayer;
+}
+
+static Int_t lbl;
+
+Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) {
+  //--------------------------------------------------------------------
+  //This functions reconstructs ITS tracks
+  //--------------------------------------------------------------------
+  TFile *in=(TFile*)inp;
+  TDirectory *savedir=gDirectory; 
+
+  if (!in->IsOpen()) {
+     cerr<<"AliITStrackerV2::Clusters2Tracks(): ";
+     cerr<<"file with TPC tracks is not open !\n";
+     return 1;
+  }
+
+  if (!out->IsOpen()) {
+     cerr<<"AliITStrackerV2::Clusters2Tracks(): ";
+     cerr<<"file for ITS tracks is not open !\n";
+     return 2;
+  }
+
+  in->cd();
+  TTree *tpcTree=(TTree*)in->Get("TreeT");
+  if (!tpcTree) {
+     cerr<<"AliITStrackerV2::Clusters2Tracks() ";
+     cerr<<"can't get a tree with TPC tracks !\n";
+     return 3;
+  }
+
+  AliTPCtrack *itrack=new AliTPCtrack; 
+  tpcTree->SetBranchAddress("tracks",&itrack);
+
+  out->cd();
+  TTree itsTree("TreeT","Tree with ITS tracks");
+  AliITStrackV2 *otrack=&fBestTrack;
+  itsTree.Branch("tracks","AliITStrackV2",&otrack,32000,0);
+
+  Int_t ntrk=0;
+
+  Int_t nentr=(Int_t)tpcTree->GetEntries();
+  for (Int_t i=0; i<nentr; i++) {
+
+    if (!tpcTree->GetEvent(i)) continue;
+    Int_t tpcLabel=itrack->GetLabel(); //save the TPC track label
+lbl=tpcLabel;
+
+#ifdef DEBUG
+if (TMath::Abs(tpcLabel)!=LAB) continue;
+cout<<tpcLabel<<" *****************\n";
+#endif
+
+    try {
+      ResetTrackToFollow(AliITStrackV2(*itrack));
+    } catch (const Char_t *msg) {
+      cerr<<msg<<endl;
+      continue;
+    }
+    ResetBestTrack();
+
+    fYV=gRandom->Gaus(0.,kSigmaYV); fZV=gRandom->Gaus(0.,kSigmaZV);
+
+    Double_t r2=fTrackToFollow.GetX()*fTrackToFollow.GetX();
+    Double_t x0=0.082/21.82*2.33*(45*45+40*40)/r2+2.000/41.28*0.03*75*75/r2;
+    fTrackToFollow.Improve(x0,fYV,fZV);
+
+    //Double_t xk=77.2;
+    Double_t xk=88.;
+    fTrackToFollow.PropagateTo(xk,0.,0.); //Ne if it's still there
+
+    xk-=0.005;
+    fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+    xk-=0.02;
+    fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);   //Kevlar
+       xk-=2.0;
+       fTrackToFollow.PropagateTo(xk, 2.0/41.28*0.029, 2.0*0.029);//Nomex
+    xk-=0.02;
+    fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);   //Kevlar
+    xk-=0.005;
+    fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+
+    xk-=14.5;
+    fTrackToFollow.PropagateTo(xk,0.,0.); //C02
+
+    xk -=0.005;
+    fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7);    //Al    
+    xk -=0.005;
+    fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71);  //Tedlar
+    xk -=0.02;
+    fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);    //Kevlar
+       xk -=0.5;
+       fTrackToFollow.PropagateTo(xk, 0.5/41.28*.029, 0.5*0.029);  //Nomex    
+    xk -=0.02;
+    fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);    //Kevlar
+    xk -=0.005;
+    fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71);  //Tedlar
+    xk -=0.005;
+    fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7);    //Al    
+
+
+    for (FollowProlongation(); fI<kMaxLayer; fI++) {
+       while (TakeNextProlongation()) FollowProlongation();
+    }
+
+#ifdef DEBUG
+cout<<fBestTrack.GetNumberOfClusters()<<" number of clusters\n\n";
+#endif
+
+    if (fBestTrack.GetNumberOfClusters() >= kMaxLayer-kLayersToSkip) {
+       cerr << ++ntrk << "                \r";
+       fBestTrack.SetLabel(tpcLabel);
+       UseClusters(&fBestTrack);
+       itsTree.Fill();
+    }
+
+  }
+
+  itsTree.Write();
+  savedir->cd();
+  cerr<<"Number of TPC tracks: "<<nentr<<endl;
+  cerr<<"Number of prolonged tracks: "<<ntrk<<endl;
+
+  delete itrack;
+
+  return 0;
+}
+
+
+AliCluster *AliITStrackerV2::GetCluster(Int_t index) const {
+  //--------------------------------------------------------------------
+  //       Return pointer to a given cluster
+  //--------------------------------------------------------------------
+  Int_t l=(index & 0xf0000000) >> 28;
+  Int_t c=(index & 0x0fffffff) >> 00;
+  return fLayers[l].GetCluster(c);
+}
+
+
+void AliITStrackerV2::FollowProlongation() {
+  //--------------------------------------------------------------------
+  //This function finds a track prolongation 
+  //--------------------------------------------------------------------
+  Int_t tryAgain=kLayersToSkip;
+
+  while (fI) {
+    fI--;
+
+#ifdef DEBUG
+cout<<fI<<' ';
+#endif
+
+    AliITSlayer &layer=fLayers[fI];
+    AliITStrackV2 &track=fTracks[fI];
+
+    if (fI==3 || fI==1) {
+      Double_t rs=0.5*(fLayers[fI+1].GetR() + layer.GetR());
+      Double_t ds=0.034; if (fI==3) ds=0.039;
+      Double_t dx0r=ds/21.82*2.33, dr=ds*2.33;
+      fTrackToFollow.Propagate(fTrackToFollow.GetAlpha(),rs,1*dx0r,dr);
+    }
+
+    //find intersection
+    Double_t x,y,z;  
+    if (!fTrackToFollow.GetGlobalXYZat(layer.GetR(),x,y,z)) {
+     cerr<<"AliITStrackerV2::FollowProlongation: failed to estimate track !\n";
+     break;
+    }
+    Double_t phi=TMath::ATan2(y,x);
+    Double_t d=layer.GetThickness(phi,z);
+    Int_t idet=layer.FindDetectorIndex(phi,z);
+    if (idet<0) {
+    cerr<<"AliITStrackerV2::FollowProlongation: failed to find a detector !\n";
+      break;
+    }
+
+    //propagate to the intersection
+    const AliITSdetector &det=layer.GetDetector(idet);
+    if (!fTrackToFollow.Propagate(det.GetPhi(),det.GetR(),1*d/21.82*2.33,d*2.33))
+    {
+      cerr<<"AliITStrackerV2::FollowProlongation: propagation failed !\n";
+      break;
+    }
+    fTrackToFollow.SetDetectorIndex(idet);
+
+    //Select possible prolongations and store the current track estimation
+    track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow);
+    Double_t dz=3*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[fI]);
+    if (dz<0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl());
+    if (dz > kMaxRoad/4) {
+      //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Z !\n";
+      break;
+    }
+    Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[fI]);
+    if (dy<0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp());
+    if (dy > kMaxRoad/4) {
+      //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Y !\n";
+      break;
+    }
+
+    Double_t zmin=track.GetZ() - dz;
+    Double_t zmax=track.GetZ() + dz;
+    Double_t ymin=track.GetY() + layer.GetR()*det.GetPhi() - dy;
+    Double_t ymax=track.GetY() + layer.GetR()*det.GetPhi() + dy;
+    if (ymax>layer.GetR()*2*TMath::Pi()) {
+       ymax-=layer.GetR()*2*TMath::Pi();
+       ymin-=layer.GetR()*2*TMath::Pi();
+    }
+    layer.SelectClusters(zmin,zmax,ymin,ymax);
+
+    //if (1/TMath::Abs(track.Get1Pt())<0.2)
+    //cout<<fI<<' '<<1/TMath::Abs(track.Get1Pt())<<' '
+    //    <<dy<<' '<<dz<<' '<<layer.InRoad()<<endl;
+
+    //take another prolongation
+    if (!TakeNextProlongation()) if (!tryAgain--) break;
+    tryAgain=kLayersToSkip;
+
+  } 
+
+  //deal with the best track
+  Int_t ncl=fTrackToFollow.GetNumberOfClusters();
+  Int_t nclb=fBestTrack.GetNumberOfClusters();
+  if (ncl)
+  if (ncl >= nclb) {
+     Double_t chi2=fTrackToFollow.GetChi2();
+     if (chi2/ncl < kChi2PerCluster) {        
+        if (ncl > nclb || chi2 < fBestTrack.GetChi2()) {
+           ResetBestTrack();
+        }
+     }
+  }
+
+  if (fI) fI++;
+}
+
+
+Int_t AliITStrackerV2::TakeNextProlongation() {
+  //--------------------------------------------------------------------
+  //This function takes another track prolongation 
+  //--------------------------------------------------------------------
+  //Double_t m[20];
+  Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
+
+  AliITSlayer &layer=fLayers[fI];
+  AliITStrackV2 &t=fTracks[fI];
+
+  Double_t dz=4*TMath::Sqrt(t.GetSigmaZ2() + kSigmaZ2[fI]);
+  Double_t dy=4*TMath::Sqrt(t.GetSigmaY2() + kSigmaY2[fI]);
+
+  const AliITSclusterV2 *c=0; Int_t ci=-1;
+  Double_t chi2=12345.;
+  while ((c=layer.GetNextCluster(ci))!=0) {
+    //if (fI==0)
+    //if (c->GetLabel(0)!=TMath::Abs(lbl)) continue; //88888888888888888888
+    Int_t idet=c->GetDetectorIndex();
+
+    if (t.GetDetectorIndex()!=idet) {
+       const AliITSdetector &det=layer.GetDetector(idet);
+       if (!t.Propagate(det.GetPhi(),det.GetR(),0.,0.)) {
+         cerr<<"AliITStrackerV2::TakeNextProlongation: propagation failed !\n";
+         continue;
+       }
+       t.SetDetectorIndex(idet);
+
+#ifdef DEBUG
+cout<<fI<<" change detector !\n";
+#endif
+
+    }
+
+    if (TMath::Abs(t.GetZ() - c->GetZ()) > dz) continue;
+    if (TMath::Abs(t.GetY() - c->GetY()) > dy) continue;
+
+    //m[0]=fYV; m[1]=fZV;
+    //chi2=t.GetPredictedChi2(c,m,d/21.82*2.33);
+    chi2=t.GetPredictedChi2(c);
+
+    if (chi2<kMaxChi2) break;
+  }
+
+#ifdef DEBUG
+cout<<fI<<" chi2="<<chi2<<' '<<t.GetY()<<' '<<t.GetZ()<<' '<<dy<<' '<<dz<<endl;
+#endif
+
+  if (chi2>=kMaxChi2) return 0;
+  if (!c) return 0;
+
+  ResetTrackToFollow(t);
+
+  //if (!fTrackToFollow.Update(m,chi2,(fI<<28)+ci)) {
+  if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) {
+     cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
+     return 0;
+  }
+  fTrackToFollow.Improve(d/21.82*2.33,fYV,fZV);
+
+#ifdef DEBUG
+cout<<"accepted lab="<<c->GetLabel(0)<<' '
+    <<fTrackToFollow.GetNumberOfClusters()<<' '
+    <<fTrackToFollow.GetY()<<' '<<fTrackToFollow.GetZ()<<endl<<endl;
+#endif
+
+  return 1;
+}
+
+
+
+
+AliITStrackerV2::AliITSlayer::AliITSlayer() {
+  //--------------------------------------------------------------------
+  //default AliITSlayer constructor
+  //--------------------------------------------------------------------
+  fN=0;
+  fDetectors=0;
+}
+
+AliITStrackerV2::AliITSlayer::
+AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd) {
+  //--------------------------------------------------------------------
+  //main AliITSlayer constructor
+  //--------------------------------------------------------------------
+  fR=r; fPhiOffset=p; fZOffset=z;
+  fNladders=nl; fNdetectors=nd;
+  fDetectors=new AliITSdetector[fNladders*fNdetectors];
+
+  fN=0;
+  fI=0;
+}
+
+AliITStrackerV2::AliITSlayer::~AliITSlayer() {
+  //--------------------------------------------------------------------
+  // AliITSlayer destructor
+  //--------------------------------------------------------------------
+  delete[] fDetectors;
+  for (Int_t i=0; i<fN; i++) delete fClusters[i];
+}
+
+Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSclusterV2 *c) {
+  //--------------------------------------------------------------------
+  //This function adds a cluster to this layer
+  //--------------------------------------------------------------------
+  if (fN==kMaxClusterPerLayer) {
+ cerr<<"AliITStrackerV2::AliITSlayer::InsertCluster(): Too many clusters !\n"; 
+   return 1;
+  }
+
+  if (fN==0) {fClusters[fN++]=c; return 0;}
+  Int_t i=FindClusterIndex(c->GetZ());
+  memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliITSclusterV2*));
+  fClusters[i]=c; fN++;
+
+  return 0;
+}
+
+Int_t AliITStrackerV2::AliITSlayer::FindClusterIndex(Double_t z) const {
+  //--------------------------------------------------------------------
+  // This function returns the index of the nearest cluster 
+  //--------------------------------------------------------------------
+  if (fN==0) return 0;
+  if (z <= fClusters[0]->GetZ()) return 0;
+  if (z > fClusters[fN-1]->GetZ()) return fN;
+  Int_t b=0, e=fN-1, m=(b+e)/2;
+  for (; b<e; m=(b+e)/2) {
+    if (z > fClusters[m]->GetZ()) b=m+1;
+    else e=m; 
+  }
+  return m;
+}
+
+void AliITStrackerV2::AliITSlayer::
+SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
+  //--------------------------------------------------------------------
+  // This function sets the "window"
+  //--------------------------------------------------------------------
+  fI=FindClusterIndex(zmin); fZmax=zmax;
+  fYmin=ymin; fYmax=ymax;
+}
+
+const AliITSclusterV2 *AliITStrackerV2::AliITSlayer::GetNextCluster(Int_t &ci){
+  //--------------------------------------------------------------------
+  // This function returns clusters within the "window" 
+  //--------------------------------------------------------------------
+  const AliITSclusterV2 *cluster=0;
+  for (Int_t i=fI; i<fN; i++) {
+    const AliITSclusterV2 *c=fClusters[i];
+    if (c->GetZ() > fZmax) break;
+    if (c->IsUsed()) continue;
+    const AliITSdetector &det=GetDetector(c->GetDetectorIndex());    
+    Double_t y=fR*det.GetPhi() + c->GetY();
+
+    if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
+    if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
+
+    if (y<fYmin) continue;
+    if (y>fYmax) continue;
+    cluster=c; ci=i;
+    fI=i+1;
+    break; 
+  }
+  return cluster;
+}
+
+Int_t AliITStrackerV2::AliITSlayer::
+FindDetectorIndex(Double_t phi, Double_t z) const {
+  //--------------------------------------------------------------------
+  //This function finds the detector crossed by the track
+  //--------------------------------------------------------------------
+  Double_t dphi=phi-fPhiOffset;
+  if      (dphi <  0) dphi += 2*TMath::Pi();
+  else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
+  Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
+
+  Double_t dz=fZOffset-z;
+  Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
+
+  if (np>=fNladders) np-=fNladders;
+  if (np<0)          np+=fNladders;
+
+#ifdef DEBUG
+cout<<np<<' '<<nz<<endl;
+#endif
+
+  return np*fNdetectors + nz;
+}
+
+Double_t 
+AliITStrackerV2::AliITSlayer::GetThickness(Double_t phi, Double_t z) const {
+  //--------------------------------------------------------------------
+  //This function returns the thickness of this layer
+  //--------------------------------------------------------------------
+  //-pi<phi<+pi
+  if (3 <fR&&fR<8 ) return 1.1*0.096;
+  if (13<fR&&fR<26) return 1.1*0.088;
+  if (37<fR&&fR<41) return 1.1*0.085;
+  return 1.1*0.081;
+}
+
+
+Double_t AliITStrackerV2::GetEffectiveThickness(Double_t phi,Double_t z) const
+{
+  //--------------------------------------------------------------------
+  //Returns the thickness between the current layer and the vertex
+  //--------------------------------------------------------------------
+  //-pi<phi<+pi
+  Double_t d=0.1;
+
+  Double_t xn=fLayers[fI].GetR();
+  for (Int_t i=0; i<fI; i++) {
+    Double_t xi=fLayers[i].GetR();
+    d+=fLayers[i].GetThickness(phi,z)*xi*xi;
+  }
+
+  if (fI>1) {
+    Double_t xi=0.5*(fLayers[1].GetR()+fLayers[2].GetR());
+    d+=0.034*xi*xi;
+  }
+
+  if (fI>3) {
+    Double_t xi=0.5*(fLayers[3].GetR()+fLayers[4].GetR());
+    d+=0.039*xi*xi;
+  }
+  return d/(xn*xn);
+}
+
+
+
+Int_t AliITStrackerV2::AliITSlayer::InRoad() const {
+  //--------------------------------------------------------------------
+  // This function returns number of clusters within the "window" 
+  //--------------------------------------------------------------------
+  Int_t ncl=0;
+  for (Int_t i=fI; i<fN; i++) {
+    const AliITSclusterV2 *c=fClusters[i];
+    if (c->GetZ() > fZmax) break;
+    //if (c->IsUsed()) continue;
+    const AliITSdetector &det=GetDetector(c->GetDetectorIndex());    
+    Double_t y=fR*det.GetPhi() + c->GetY();
+
+    if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
+    if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
+
+    if (y<fYmin) continue;
+    if (y>fYmax) continue;
+    ncl++;
+  }
+  return ncl;
+}
+
+
+
diff --git a/ITS/AliITStrackerV2.h b/ITS/AliITStrackerV2.h
new file mode 100644 (file)
index 0000000..d01c771
--- /dev/null
@@ -0,0 +1,111 @@
+#ifndef ALIITSTRACKER_H
+#define ALIITSTRACKER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                          ITS tracker
+//
+//       Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
+//-------------------------------------------------------------------------
+#include "AliTracker.h"
+#include "AliITSrecoV2.h"
+#include "AliITStrackV2.h"
+
+class AliITSclusterV2;
+class AliITSgeom;
+class TFile;
+
+
+//-------------------------------------------------------------------------
+class AliITStrackerV2 : public AliTracker {
+public:
+  AliITStrackerV2():AliTracker(){}
+  AliITStrackerV2(const AliITSgeom *geom) throw (const Char_t *);
+
+  AliCluster *GetCluster(Int_t index) const;
+  Int_t Clusters2Tracks(const TFile *in, TFile *out);
+  Int_t PropagateBack(const TFile *in, TFile *out) {return 0;}
+
+private:
+
+  Double_t GetEffectiveThickness(Double_t phi, Double_t z) const;
+
+  void  FollowProlongation();
+  Int_t TakeNextProlongation();
+
+  void ResetBestTrack() {
+     fBestTrack.~AliITStrackV2();
+     new(&fBestTrack) AliITStrackV2(fTrackToFollow);
+  }
+
+  void ResetTrackToFollow(const AliITStrackV2 &t) {
+     fTrackToFollow.~AliITStrackV2();
+     new(&fTrackToFollow) AliITStrackV2(t);
+  }
+
+class AliITSdetector {
+private:
+  Double_t fR;    // polar coordinates 
+  Double_t fPhi;  // of this detector
+
+public:
+  AliITSdetector(){}
+  AliITSdetector(Double_t r,Double_t phi) {fR=r; fPhi=phi;}
+
+  void *operator new(size_t s,AliITSdetector *p) {return p;}
+
+  Double_t GetR()   const {return fR;}
+  Double_t GetPhi() const {return fPhi;}
+};
+
+class AliITSlayer {
+  Double_t fR;                // mean radius of this layer
+  Double_t fPhiOffset;        // offset of the first detector in Phi
+  Int_t fNladders;            // number of ladders
+  Double_t fZOffset;          // offset of the first detector in Z
+  Int_t fNdetectors;          // detectors/ladder
+  AliITSdetector *fDetectors; // array of detectors
+
+  Int_t fN;                   // number of clusters
+  AliITSclusterV2 *fClusters[kMaxClusterPerLayer];    // pointers to clusters
+
+  Double_t fZmax;      //       edges
+  Double_t fYmin;      //      of  the
+  Double_t fYmax;      //      "window"
+  Int_t fI;            // index of the current cluster within the "window"
+
+  Int_t FindClusterIndex(Double_t z) const;
+
+public:
+  AliITSlayer();
+  AliITSlayer(Double_t r, Double_t p, Double_t z, Int_t nl, Int_t nd);
+ ~AliITSlayer();
+  Int_t InsertCluster(AliITSclusterV2 *c);
+  void SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin,Double_t ymax);
+  const AliITSclusterV2 *GetNextCluster(Int_t &ci);
+
+  void *operator new(size_t s, AliITSlayer *p) {return p;}
+
+  Double_t GetR() const {return fR;}
+  AliITSclusterV2 *GetCluster(Int_t i) const {return fClusters[i];} 
+  AliITSdetector &GetDetector(Int_t n) const { return fDetectors[n]; }
+  Int_t FindDetectorIndex(Double_t phi, Double_t z) const;
+  Double_t GetThickness(Double_t phi, Double_t z) const;
+
+  Int_t InRoad() const ;
+};
+
+  Int_t fI;                       // index of the current layer
+  static AliITSlayer fLayers[kMaxLayer]; // ITS layers
+  AliITStrackV2 fTracks[kMaxLayer]; // track estimations at the ITS layers
+  AliITStrackV2 fBestTrack;         // "best" track 
+  AliITStrackV2 fTrackToFollow;     // followed track
+
+  Double_t fYV;                   // Y-coordinate of the primary vertex
+  Double_t fZV;                   // Z-coordinate of the primary vertex
+};
+
+
+
+#endif
index 6039220..c7dbe3a 100644 (file)
 // This class will always be for ITS only
 #pragma link C++ class  AliITSvtest+;
 
+#pragma link C++ class AliITSclusterV2+;
+#pragma link C++ class AliITStrackV2+;
+#pragma link C++ class AliITStrackerV2+;
+
 #endif
index 68f83df..d8a2b52 100644 (file)
@@ -39,7 +39,8 @@ SRCS          = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \
                AliITSstatistics.cxx AliITSstatistics2.cxx \
                 AliITStrack.cxx AliITStracking.cxx AliITSiotrack.cxx \
                AliITSRad.cxx \
-               AliITSvtest.cxx         
+               AliITSvtest.cxx \
+                AliITSclusterV2.cxx AliITStrackV2.cxx AliITStrackerV2.cxx      
 #              AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \
 #              AliITStrack.cxx