--- /dev/null
+#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;
+}
+
+
--- /dev/null
+#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;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+#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;
+}
--- /dev/null
+/**************************************************************************
+ * 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)
--- /dev/null
+#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
+
+
--- /dev/null
+#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
--- /dev/null
+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;
+}
--- /dev/null
+/**************************************************************************
+ * 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;
+}
+*/
+
--- /dev/null
+#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
+
+
--- /dev/null
+/**************************************************************************
+ * 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;
+}
+
+
+
--- /dev/null
+#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
// 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
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
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff