Updated V2 stream of tracking (Yu.Belikov). The new long waited features are: 1)...
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 8 Nov 2001 16:36:33 +0000 (16:36 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 8 Nov 2001 16:36:33 +0000 (16:36 +0000)
26 files changed:
ITS/AliITSComparisonV2.C
ITS/AliITSFindClustersV2.C
ITS/AliITSFindTracksV2.C
ITS/AliITSrecoV2.h
ITS/AliITStrackV2.cxx
ITS/AliITStrackV2.h
ITS/AliITStrackerV2.cxx
ITS/AliITStrackerV2.h
ITS/AliV0Comparison.C [new file with mode: 0644]
ITS/AliV0FindVertices.C [new file with mode: 0644]
ITS/AliV0vertex.cxx [new file with mode: 0644]
ITS/AliV0vertex.h [new file with mode: 0644]
ITS/AliV0vertexer.cxx [new file with mode: 0644]
ITS/AliV0vertexer.h [new file with mode: 0644]
ITS/ITSLinkDef.h
ITS/Makefile
ITS/libITS.pkg
STEER/AliKalmanTrack.h
STEER/AliTracker.cxx
STEER/AliTracker.h
TPC/AliTPCComparison.C
TPC/AliTPCFindTracks.C
TPC/AliTPCtrack.cxx
TPC/AliTPCtrack.h
TPC/AliTPCtracker.cxx
TPC/AliTPCtracker.h

index ef036cf..2f82cc2 100644 (file)
 #endif
 
 struct GoodTrackITS {
-  Int_t event;
   Int_t lab;
   Int_t code;
   Float_t px,py,pz;
   Float_t x,y,z;
 };
 
-Int_t AliITSComparisonV2() {
-   Int_t good_tracks_its(GoodTrackITS *gt, Int_t max);
-
+Int_t AliITSComparisonV2(Int_t event=0) {
    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_ITS_0");
-   if (!tracktree) {cerr<<"Can't get a tree with ITS tracks !\n"; return 4;}
-   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);
+   const Int_t MAX=15000;
+   Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const Int_t event);
+
+   Int_t nentr=0; TObjArray tarray(2000);
+   {/* Load tracks */ 
+     TFile *tf=TFile::Open("AliITStracksV2.root");
+     if (!tf->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 3;}
+     char tname[100]; sprintf(tname,"TreeT_ITS_%d",event);
+     TTree *tracktree=(TTree*)tf->Get(tname);
+     if (!tracktree) {cerr<<"Can't get a tree with ITS tracks !\n"; return 4;}
+     TBranch *tbranch=tracktree->GetBranch("tracks");
+     nentr=(Int_t)tracktree->GetEntries();
+     for (Int_t i=0; i<nentr; i++) {
+        AliITStrackV2 *iotrack=new AliITStrackV2;
+        tbranch->SetAddress(&iotrack);
+        tracktree->GetEvent(i);
+        tarray.AddLast(iotrack);
+        /*if (itsLabel != 1234) continue;
+        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;
+       }*/
+     }
+     delete tracktree; //Thanks to Mariana Bondila
+     tf->Close();
    }
-   delete tracktree; //Thanks to Mariana Bondila
-   tf->Close();
-   delete geom; //Thanks to Mariana Bondila
-   cf->Close();
 
-/////////////////////////////////////////////////////////////////////////
-   const Int_t MAX=15000;
+   /* Generate a list of "good" tracks */
    GoodTrackITS gt[MAX];
    Int_t ngood=0;
    ifstream in("good_tracks_its");
    if (in) {
       cerr<<"Reading good tracks...\n";
-      while (in>>gt[ngood].event>>gt[ngood].lab>>gt[ngood].code>>
+      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++;
@@ -94,11 +78,11 @@ Int_t AliITSComparisonV2() {
       if (!in.eof()) cerr<<"Read error (good_tracks_its) !\n";
    } else {
       cerr<<"Marking good tracks (this will take a while)...\n";
-      ngood=good_tracks_its(gt,MAX);
+      ngood=good_tracks_its(gt,MAX,event);
       ofstream out("good_tracks_its");
       if (out) {
        for (Int_t ngd=0; ngd<ngood; ngd++)
-         out<<gt[ngd].event<<' '<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
+         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";
@@ -246,7 +230,7 @@ Int_t AliITSComparisonV2() {
    return 0;
 }
 
-Int_t good_tracks_its(GoodTrackITS *gt, Int_t max) {
+Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const Int_t event) {
    if (gAlice) {delete gAlice; gAlice=0;}
 
    TFile *file=TFile::Open("galice.root");
@@ -256,7 +240,7 @@ Int_t good_tracks_its(GoodTrackITS *gt, Int_t max) {
      exit(5);
    }
 
-   Int_t np=gAlice->GetEvent(0);
+   Int_t np=gAlice->GetEvent(event);
 
    Int_t *good=new Int_t[np];
    Int_t k;
@@ -275,7 +259,8 @@ Int_t good_tracks_its(GoodTrackITS *gt, Int_t max) {
    if (!cf->IsOpen()){
       cerr<<"Can't open AliITSclustersV2.root !\n"; exit(6);
    }
-   TTree *cTree=(TTree*)cf->Get("TreeC_ITS_0");
+   char cname[100]; sprintf(cname,"TreeC_ITS_%d",event);
+   TTree *cTree=(TTree*)cf->Get(cname);
    if (!cTree) {
       cerr<<"Can't get cTree !\n"; exit(7);
    }
@@ -315,11 +300,10 @@ Int_t good_tracks_its(GoodTrackITS *gt, Int_t max) {
    }
    Int_t nt=0;
    Double_t px,py,pz,x,y,z;
-   Int_t code,lab,event;
-   while (in>>event>>lab>>code>>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].event=event;
       gt[nt].lab=lab;
       gt[nt].code=p->GetPdgCode();
 //**** px py pz - in global coordinate system
index bc82fa9..6f9d151 100644 (file)
@@ -82,7 +82,6 @@ Int_t AliITSFindClustersV2() {
    geom->Write();
 
    TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
-   //TTree *cTree=new TTree("cTree","ITS clusters");
    TTree *cTree=new TTree("TreeC_ITS_0","ITS clusters");
    cTree->Branch("Clusters",&clusters);
 
@@ -105,6 +104,8 @@ Int_t AliITSFindClustersV2() {
 
    cerr<<"Number of entries: "<<nentr<<endl;
 
+   Float_t lp[5]; Int_t lab[6]; //Why can't it be inside a loop ?
+
    for (Int_t i=0; i<nentr; i++) {
        points->Clear();
        pTree->GetEvent(i);
@@ -117,13 +118,13 @@ Int_t AliITSFindClustersV2() {
        nclusters+=ncl;
        for (Int_t j=0; j<ncl; j++) {
           AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
-          Float_t lp[5];
+          //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]; 
+          //Int_t lab[6]; 
           lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
           lab[3]=ndet;
 
@@ -143,7 +144,6 @@ Int_t AliITSFindClustersV2() {
           new(cl[j]) AliITSclusterV2(lab,lp);
        }
        cTree->Fill(); clusters->Delete();
-       points->Delete();
    }
    cTree->Write();
 
index aa3ef1d..00ac80a 100644 (file)
@@ -1,5 +1,6 @@
 #ifndef __CINT__
   #include <iostream.h>
+  #include "AliITSgeom.h"
   #include "AliITStrackerV2.h"
 
   #include "TFile.h"
@@ -22,6 +23,12 @@ Int_t AliITSFindTracksV2() {
 
    TStopwatch timer;
    AliITStrackerV2 tracker(geom);
+
+   //Double_t xyz[]={0.,0.,0.}; tracker.SetVertex(xyz);  //primary vertex
+   //Int_t flag[]={1};                                   //some default flags
+   //flag[0]= 0; tracker.SetupFirstPass(flag);           //no constraint
+   //flag[0]=-1; tracker.SetupSecondPass(flag);          //skip second pass
+
    Int_t rc=tracker.Clusters2Tracks(in,out);
    timer.Stop(); timer.Print();
 
index 93603a7..ef11775 100644 (file)
@@ -24,9 +24,9 @@
      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 kChi2PerCluster=7.;//5.;
+   const Double_t kMaxChi2=15.;//12;
+   const Double_t kMaxRoad=3.0;
 
    const Double_t kSigmaYV=0.005e+0;
    const Double_t kSigmaZV=0.010e+0;
index aa7ad01..3dfb99f 100644 (file)
@@ -30,7 +30,7 @@
 
 ClassImp(AliITStrackV2)
 
-const Int_t kWARN=1;
+const Int_t kWARN=5;
 
 //____________________________________________________________________________
 AliITStrackV2::AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *) {
@@ -42,7 +42,9 @@ AliITStrackV2::AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *) {
   SetNumberOfClusters(0);
   //SetConvConst(t.GetConvConst());
 
-  fdEdx  = 0.;
+  fdEdx  = t.GetdEdx();
+  SetMass(t.GetMass());
+
   fAlpha = t.GetAlpha();
   if      (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
   else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
@@ -88,19 +90,35 @@ AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : AliKalmanTrack(t) {
   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;
+  AliITStrackV2 *t=(AliITStrackV2*)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;
 }
+*/
+//_____________________________________________________________________________
+Int_t AliITStrackV2::Compare(const TObject *o) const {
+  //-----------------------------------------------------------------
+  // This function compares tracks according to the their curvature
+  //-----------------------------------------------------------------
+  AliITStrackV2 *t=(AliITStrackV2*)o;
+
+  Double_t p2=1./(Get1Pt()*Get1Pt());
+  Double_t b2=p2/(p2 + GetMass()*GetMass());
+  Double_t po2=1./(t->Get1Pt()*t->Get1Pt());
+  Double_t bo2=po2/(po2 + t->GetMass()*t->GetMass());
+  if (p2*b2>po2*bo2) return -1;
+  else if (p2*b2<po2*bo2) return 1;
+  return 0;
+}
 
 //_____________________________________________________________________________
 void AliITStrackV2::GetExternalCovariance(Double_t cc[15]) const {
@@ -118,12 +136,12 @@ void AliITStrackV2::GetExternalCovariance(Double_t cc[15]) const {
 }
 
 //____________________________________________________________________________
-Int_t AliITStrackV2::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm) {
+Int_t AliITStrackV2::PropagateToVertex(Double_t x0,Double_t rho) {
   //------------------------------------------------------------------
   //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);   
+  Propagate(fAlpha,xv,0.,0.);   
   return 0;
 }
 
@@ -136,8 +154,9 @@ GetGlobalXYZat(Double_t xk, Double_t &x, Double_t &y, Double_t &z) const {
   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";
+    Int_t n=GetNumberOfClusters();
+    if (n>kWARN) 
+      cerr<<n<<" AliITStrackV2::GetGlobalXYZat: Propagation failed !\n";
     return 0;
   }
 
@@ -164,9 +183,10 @@ Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const
   r00+=fC00; r01+=fC10; r11+=fC11;
 
   Double_t det=r00*r11 - r01*r01;
-  if (TMath::Abs(det) < 1.e-10) {
+  if (TMath::Abs(det) < 1.e-30) {
     Int_t n=GetNumberOfClusters();
-    if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
+    if (n>kWARN) 
+       cerr<<n<<" AliKalmanTrack::GetPredictedChi2: Singular matrix !\n";
     return 1e10;
   }
   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
@@ -178,7 +198,7 @@ Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const
 
 //_____________________________________________________________________________
 Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c,Double_t *m,
-Double_t x0, Double_t pm) const {
+Double_t x0) const {
   //-----------------------------------------------------------------
   // This function calculates a chi2 increment with a vertex contraint 
   //-----------------------------------------------------------------
@@ -209,7 +229,7 @@ Double_t x0, Double_t pm) const {
   //x0=0.;
   if (x0!=0.) {
      Double_t pp2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
-     Double_t beta2=pp2/(pp2 + pm*pm);
+     Double_t beta2=pp2/(pp2 + GetMass()*GetMass());
      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());
@@ -235,9 +255,10 @@ Double_t x0, Double_t pm) const {
   TMatrixD R(tmp,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed,H)); R+=V;
   
   Double_t det=R.Determinant();
-  if (TMath::Abs(det) < 1.e-25) {
+  if (TMath::Abs(det) < 1.e-30) {
     Int_t n=GetNumberOfClusters();
-    if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Singular matrix !\n";
+    if (n>kWARN) 
+       cerr<<n<<" AliITStrackV2::GetPredictedChi2: Singular matrix !\n";
     return 1e10;
   }
 
@@ -250,7 +271,7 @@ Double_t x0, Double_t pm) const {
 
 //____________________________________________________________________________
 Int_t 
-AliITStrackV2::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) {
+AliITStrackV2::PropagateTo(Double_t xk,Double_t x0,Double_t rho) {
   //------------------------------------------------------------------
   //This function propagates a track
   //------------------------------------------------------------------
@@ -258,7 +279,8 @@ AliITStrackV2::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) {
   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";
+    if (n>kWARN) 
+       cerr<<n<<" AliITStrackV2::PropagateTo: Propagation failed !\n";
     return 0;
   }
 
@@ -311,7 +333,7 @@ AliITStrackV2::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) {
   fX=x2;
 
   Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
-  Double_t beta2=p2/(p2 + pm*pm);
+  Double_t beta2=p2/(p2 + GetMass()*GetMass());
 
   //Multiple scattering******************
   //x0=0.;
@@ -329,7 +351,7 @@ AliITStrackV2::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) {
      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);
+     fP4*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
   }
 
   if (!Invariant()) {cout<<"Propagate !\n"; return 0;}
@@ -363,13 +385,7 @@ Int_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, UInt_t index) {
 
   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;
@@ -492,15 +508,19 @@ Int_t AliITStrackV2::Invariant() const {
   //------------------------------------------------------------------
   // This function is for debugging purpose only
   //------------------------------------------------------------------
+  Int_t n=GetNumberOfClusters();
+  
   //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;}
+  //if (fP1*fP4<0) {
+  //   if (n>kWARN) cout<<"fP1*fP4="<<fP1*fP4<<' '<<fP1<<endl; return 0;}
+  
+  if (TMath::Abs(fP2)>=1) {if (n>kWARN) cout<<"fP2="<<fP2<<endl; return 0;}
+
+  if (fC00<=0) {if (n>kWARN) cout<<"fC00="<<fC00<<endl; return 0;}
+  if (fC11<=0) {if (n>kWARN) cout<<"fC11="<<fC11<<endl; return 0;}
+  if (fC22<=0) {if (n>kWARN) cout<<"fC22="<<fC22<<endl; return 0;}
+  if (fC33<=0) {if (n>kWARN) cout<<"fC33="<<fC33<<endl; return 0;}
+  if (fC44<=0) {if (n>kWARN) cout<<"fC44="<<fC44<<endl; return 0;}
   /*
   TMatrixD m(5,5);
   m(0,0)=fC00; 
@@ -517,8 +537,7 @@ Int_t AliITStrackV2::Invariant() const {
   Double_t det=m.Determinant(); 
 
   if (det <= 0) {
-      cout<<" bad determinant "<<det<<endl;
-      m.Print(); 
+      if (n>kWARN) { cout<<" bad determinant "<<det<<endl; m.Print(); } 
       return 0;
   }
   */
@@ -526,8 +545,8 @@ Int_t AliITStrackV2::Invariant() const {
 }
 
 //____________________________________________________________________________
-Int_t AliITStrackV2::Propagate(Double_t alp, Double_t xk,
-Double_t x0,Double_t rho,Double_t pm) {
+Int_t 
+AliITStrackV2::Propagate(Double_t alp,Double_t xk,Double_t x0,Double_t rho) {
   //------------------------------------------------------------------
   //This function propagates a track
   //------------------------------------------------------------------
@@ -547,7 +566,8 @@ Double_t x0,Double_t rho,Double_t pm) {
   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";
+     if (n>kWARN) 
+        cerr<<n<<" AliITStrackV2::Propagate: Rotation failed !\n";
      return 0;
   }
 
@@ -569,7 +589,8 @@ Double_t x0,Double_t rho,Double_t pm) {
   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";
+    if (n>kWARN) 
+       cerr<<n<<" AliITStrackV2::Propagate: Propagation failed !\n";
     return 0;
   }
 
@@ -646,7 +667,7 @@ Double_t x0,Double_t rho,Double_t pm) {
   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);
+  Double_t beta2=pp2/(pp2 + GetMass()*GetMass());
 
   //Multiple scattering******************
   //x0=0.;
@@ -664,7 +685,7 @@ Double_t x0,Double_t rho,Double_t pm) {
      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);
+     fP4*=(1.- sqrt(pp2+GetMass()*GetMass())/pp2*dE);
   }
 
   if (!Invariant()) {
@@ -700,7 +721,7 @@ Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
   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);
+  Double_t beta2=p2/(p2 + GetMass()*GetMass());
   x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
 
@@ -742,7 +763,7 @@ Int_t AliITStrackV2::Improve(Double_t x0,Double_t xv,Double_t yv) {
 
   Double_t r2=fX*fX+fP0*fP0;
   Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
-  Double_t beta2=p2/(p2 + 0.14*0.14);
+  Double_t beta2=p2/(p2 + GetMass()*GetMass());
   x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
 
index c3660eb..aab0a25 100644 (file)
@@ -36,10 +36,17 @@ 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;}
+  Int_t PropagateToVertex(Double_t x0=36.66,Double_t rho=1.2e-3);
+  Int_t Propagate(Double_t alpha, Double_t xr, Double_t x0, Double_t rho);
+  Int_t PropagateTo(Double_t xr,Double_t x0=21.82,Double_t rho=2.33);
+  Int_t Update(const Double_t *m, Double_t chi2, UInt_t i);
+  Int_t Update(const AliCluster* cl,Double_t chi2,UInt_t i);
+  Int_t Improve(Double_t x0,Double_t yv,Double_t zv);
+  void SetdEdx(Double_t dedx) {fdEdx=dedx;}
+  void CookdEdx(Double_t low=0., Double_t up=1.) {}
   void SetDetectorIndex(Int_t i) {SetLabel(i);}
+  void ResetCovariance();
+  void ResetClusters() { SetChi2(0.); SetNumberOfClusters(0); }
   
   void *operator new(size_t s,void *p) { return p; }
   void *operator new(size_t s) { return ::operator new(s); }
@@ -47,7 +54,7 @@ public:
   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 GetdEdx() const {return fdEdx;}
   Double_t GetY()    const {return fP0;}
   Double_t GetZ()    const {return fP1;}
   Double_t GetSnp()  const {return fP2;}
@@ -56,33 +63,15 @@ public:
   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);
-  void ResetCovariance();
-  void ResetClusters() { SetChi2(0.); SetNumberOfClusters(0); }
-
+  Double_t 
+  GetPredictedChi2(const AliCluster *cluster, Double_t *m, Double_t x0) const;
   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)
index 1a49986..5564549 100644 (file)
 //#define DEBUG
 
 #ifdef DEBUG
-Int_t LAB=236;
+Int_t LAB=7;
 #endif
 
-extern TRandom *gRandom;
-
 AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers
 
 AliITStrackerV2::
-AliITStrackerV2(const AliITSgeom *geom, Int_t eventn) throw (const Char_t *) {
+AliITStrackerV2(const AliITSgeom *geom, Int_t eventn) throw (const Char_t *) :
+AliTracker() {
   //--------------------------------------------------------------------
   //This is the AliITStracker constructor
   //It also reads clusters from a root file
   //--------------------------------------------------------------------
-  fEventN = eventn;  //MI change add event number  - used to generate identifier 
-  fYV=fZV=0.;
+  fEventN=eventn;  //MI change add event number  - used to generate identifier 
 
   AliITSgeom *g=(AliITSgeom*)geom;
 
@@ -122,10 +119,7 @@ if (c->GetLabel(0)!=LAB) continue;
 cout<<lay-1<<' '<<lad-1<<' '<<det-1<<' '<<c->GetY()<<' '<<c->GetZ()<<endl;
 #endif
 
-         AliITSdetector &det=fLayers[lay-1].GetDetector(c->GetDetectorIndex());
-         Double_t r=det.GetR(); r=TMath::Sqrt(r*r + c->GetY()*c->GetY());
-         if (r > TMath::Abs(c->GetZ())-0.5)
-            fLayers[lay-1].InsertCluster(new AliITSclusterV2(*c));
+         fLayers[lay-1].InsertCluster(new AliITSclusterV2(*c));
        }
        clusters->Delete();
      }
@@ -137,6 +131,9 @@ cout<<lay-1<<' '<<lad-1<<' '<<det-1<<' '<<c->GetY()<<' '<<c->GetZ()<<endl;
   }
 
   fI=kMaxLayer;
+
+  fPass=0;
+  fConstraint[0]=1; fConstraint[1]=0;
 }
 
 #ifdef DEBUG
@@ -162,112 +159,121 @@ Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) {
      return 2;
   }
 
-  in->cd(); 
-  //
+  in->cd();
   char   tname[100];
-  sprintf(tname,"TreeT_TPC_%d",fEventN);
-  TTree *tpcTree=(TTree*)in->Get(tname);
-
-  if (!tpcTree) {
-     cerr<<"AliITStrackerV2::Clusters2Tracks() ";
-     cerr<<"can't get a tree with TPC tracks !\n";
-     return 3;
+  Int_t nentr=0; TObjArray itsTracks(10000);
+
+  {/* Read TPC tracks */ 
+    sprintf(tname,"TreeT_TPC_%d",fEventN);
+    TTree *tpcTree=(TTree*)in->Get(tname);
+    if (!tpcTree) {
+       cerr<<"AliITStrackerV2::Clusters2Tracks(): "
+             "can't get a tree with TPC tracks !\n";
+       return 3;
+    }
+    AliTPCtrack *itrack=new AliTPCtrack; 
+    tpcTree->SetBranchAddress("tracks",&itrack);
+    nentr=(Int_t)tpcTree->GetEntries();
+    for (Int_t i=0; i<nentr; i++) {
+       tpcTree->GetEvent(i);
+       itsTracks.AddLast(new AliITStrackV2(*itrack));
+    }
+    delete tpcTree; //Thanks to Mariana Bondila
+    delete itrack;
   }
-
-  AliTPCtrack *itrack=new AliTPCtrack; 
-  tpcTree->SetBranchAddress("tracks",&itrack);
+  itsTracks.Sort();
 
   out->cd();
-  TTree itsTree("ITSf","Tree with ITS tracks");
+
+  sprintf(tname,"TreeT_ITS_%d",fEventN);
+  TTree itsTree(tname,"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
-
+  for (fPass=0; fPass<2; fPass++) {
+     Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
+     for (Int_t i=0; i<nentr; i++) {
+       if (i%10==0) cerr<<nentr-i<<" \r";
+       AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i);
+       if (t==0) continue;           //this track has been already tracked
+       Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 #ifdef DEBUG
 lbl=tpcLabel;
 if (TMath::Abs(tpcLabel)!=LAB) continue;
 cout<<tpcLabel<<" *****************\n";
 #endif
+       try {
+         ResetTrackToFollow(*t);
+       } catch (const Char_t *msg) {
+         cerr<<msg<<endl;
+         continue;
+       }
+       ResetBestTrack();
 
-    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=80.;
-    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=61.;
-    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();
-    }
+       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;
+       if (constraint) fTrackToFollow.Improve(x0,GetY(),GetZ());
+
+       Double_t xk=80.;
+       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=61.;
+       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();
-    }
+       if (fBestTrack.GetNumberOfClusters() >= kMaxLayer-kLayersToSkip) {
+          fBestTrack.SetLabel(tpcLabel);
+         fBestTrack.CookdEdx();
+          CookLabel(&fBestTrack,0.); //For comparison only
+          itsTree.Fill();
+          UseClusters(&fBestTrack);
+          delete itsTracks.RemoveAt(i);
+       }
 
+     }
   }
-  sprintf(tname,"TreeT_ITS_%d",fEventN);
-  itsTree.Write(tname);
-  savedir->cd();
-  cerr<<"Number of TPC tracks: "<<nentr<<endl;
-  cerr<<"Number of prolonged tracks: "<<ntrk<<endl;
 
-  delete tpcTree; //Thanks to Mariana Bondila
+  itsTree.Write();
 
-  delete itrack;
+  itsTracks.Delete();
+  savedir->cd();
+  cerr<<"Number of TPC tracks: "<<nentr<<endl;
+  cerr<<"Number of prolonged tracks: "<<itsTree.GetEntries()<<endl;
 
   return 0;
 }
@@ -488,48 +494,50 @@ cout<<fI<<' ';
       Double_t rs=0.5*(fLayers[fI+1].GetR() + r);
       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);
+      fTrackToFollow.Propagate(fTrackToFollow.GetAlpha(),rs,dx0r,dr);
     }
 
     //find intersection
     Double_t x,y,z;  
     if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) {
-     cerr<<"AliITStrackerV2::FollowProlongation: failed to estimate track !\n";
+      //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";
+      //cerr<<"AliITStrackerV2::FollowProlongation: "
+      //"failed to find a detector !\n";
       break;
     }
 
     //propagate to the intersection
     const AliITSdetector &det=layer.GetDetector(idet);
     phi=det.GetPhi();
-    if (!fTrackToFollow.Propagate(phi,det.GetR(),1*d/21.82*2.33,d*2.33)) {
-      cerr<<"AliITStrackerV2::FollowProlongation: propagation failed !\n";
+    if (!fTrackToFollow.Propagate(phi,det.GetR(),d/21.82*2.33,d*2.33)) {
+      //cerr<<"AliITStrackerV2::FollowProlongation: propagation failed !\n";
       break;
     }
+    if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+1) 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) {
+    if (dz > kMaxRoad) {
       //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) {
+    if (dy > kMaxRoad) {
       //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Y !\n";
       break;
     }
-
-    Double_t zmin=track.GetZ() - dz;
+    Double_t zmin=track.GetZ() - dz; 
     Double_t zmax=track.GetZ() + dz;
     Double_t ymin=track.GetY() + r*phi - dy;
     Double_t ymax=track.GetY() + r*phi + dy;
@@ -567,6 +575,7 @@ Int_t AliITStrackerV2::TakeNextProlongation() {
 
   AliITSlayer &layer=fLayers[fI];
   AliITStrackV2 &t=fTracks[fI];
+  Int_t &constraint=fConstraint[fPass];
 
   Double_t dz=4*TMath::Sqrt(t.GetSigmaZ2() + kSigmaZ2[fI]);
   Double_t dy=4*TMath::Sqrt(t.GetSigmaY2() + kSigmaY2[fI]);
@@ -585,7 +594,8 @@ Int_t AliITStrackerV2::TakeNextProlongation() {
     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";
+         //cerr<<"AliITStrackerV2::TakeNextProlongation: "
+         //"propagation failed !\n";
          continue;
        }
        t.SetDetectorIndex(idet);
@@ -617,10 +627,10 @@ cout<<fI<<" chi2="<<chi2<<' '<<t.GetY()<<' '<<t.GetZ()<<' '<<dy<<' '<<dz<<endl;
 
   //if (!fTrackToFollow.Update(m,chi2,(fI<<28)+ci)) {
   if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) {
-     cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
+     //cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
      return 0;
   }
-  fTrackToFollow.Improve(d/21.82*2.33,fYV,fZV);
+  if (constraint) fTrackToFollow.Improve(d/21.82*2.33,GetY(),GetZ());
 
 #ifdef DEBUG
 cout<<"accepted lab="<<c->GetLabel(0)<<' '
@@ -640,7 +650,6 @@ AliITStrackerV2::AliITSlayer::AliITSlayer() {
   //--------------------------------------------------------------------
   fN=0;
   fDetectors=0;
-  for (Int_t i=0; i<kMaxClusterPerLayer; i++) fClusters[i]=0;
 }
 
 AliITStrackerV2::AliITSlayer::
@@ -651,7 +660,6 @@ AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd) {
   fR=r; fPhiOffset=p; fZOffset=z;
   fNladders=nl; fNdetectors=nd;
   fDetectors=new AliITSdetector[fNladders*fNdetectors];
-  for (Int_t i=0; i<kMaxClusterPerLayer; i++) fClusters[i]=0;
 
   fN=0;
   fI=0;
@@ -729,6 +737,7 @@ const AliITSclusterV2 *AliITStrackerV2::AliITSlayer::GetNextCluster(Int_t &ci){
     fI=i+1;
     break; 
   }
+
   return cluster;
 }
 
@@ -761,7 +770,7 @@ 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 2.2*0.096;
+  if (3 <fR&&fR<8 ) return 2.5*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;
index a8e58cf..dc7f52e 100644 (file)
@@ -21,95 +21,113 @@ class TFile;
 class AliITStrackerV2 : public AliTracker {
 public:
   AliITStrackerV2():AliTracker(){}
-  AliITStrackerV2(const AliITSgeom *geom, Int_t event=0) throw (const Char_t *);
-
+  AliITStrackerV2(const AliITSgeom *geom,Int_t event=0) 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);
+  void SetupFirstPass(Int_t *flags, Double_t *cuts=0);
+  void SetupSecondPass(Int_t *flags, Double_t *cuts=0);
+
+  class AliITSdetector {
+  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;}
+  private:
+    Double_t fR;    // polar coordinates 
+    Double_t fPhi;  // of this detector
+  };
+
+  class AliITSlayer {
+  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 zmi,Double_t zma,Double_t ymi,Double_t yma);
+    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 GetNumberOfClusters() const {return fN;}
+  private:
+    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;
+  };
 
 private:
-  Int_t fEventN;  //event number
+  void CookLabel(AliKalmanTrack *t,Float_t wrong) const;
   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);
   }
-
-  //The two subclasses are public in order to access one from the another
- public: 
-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 GetNumberOfClusters() const {return fN;}
-};
-
- private:
-  Int_t fI;                       // index of the current layer
+  Int_t fEventN;                         //event number
+  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
+  AliITStrackV2 fTracks[kMaxLayer];      // track estimations at the ITS layers
+  AliITStrackV2 fBestTrack;              // "best" track 
+  AliITStrackV2 fTrackToFollow;          // followed track
+  Int_t fPass;                           // current pass through the data 
+  Int_t fConstraint[2];                  // constraint flags
 };
 
 
+inline void AliITStrackerV2::SetupFirstPass(Int_t *flags, Double_t *cuts) {
+  // This function sets up flags and cuts for the first tracking pass   
+  //
+  //   flags[0] - vertex constaint flag                                
+  //              negative means "skip the pass"                        
+  //              0        means "no constraint"                        
+  //              positive means "normal constraint"                    
+
+   fConstraint[0]=flags[0];
+   if (cuts==0) return;
+}
+
+inline void AliITStrackerV2::SetupSecondPass(Int_t *flags, Double_t *cuts) {
+  // This function sets up flags and cuts for the second tracking pass   
+  //
+  //   flags[0] - vertex constaint flag                                
+  //              negative means "skip the pass"                        
+  //              0        means "no constraint"                        
+  //              positive means "normal constraint"                    
+
+   fConstraint[1]=flags[0];
+   if (cuts==0) return;
+}
+
+inline void AliITStrackerV2::CookLabel(AliKalmanTrack *t,Float_t wrong) const {
+  //--------------------------------------------------------------------
+  //This function "cooks" a track label. If label<0, this track is fake.
+  //--------------------------------------------------------------------
+   Int_t tpcLabel=t->GetLabel();
+   if (tpcLabel<0) return;
+   AliTracker::CookLabel(t,wrong);
+   if (tpcLabel != t->GetLabel()) t->SetLabel(-tpcLabel); 
+}
 
 #endif
diff --git a/ITS/AliV0Comparison.C b/ITS/AliV0Comparison.C
new file mode 100644 (file)
index 0000000..53bc2f3
--- /dev/null
@@ -0,0 +1,336 @@
+#ifndef __CINT__
+  #include <iostream.h>
+  #include <fstream.h>
+
+  #include "TH1.h"
+  #include "TFile.h"
+  #include "TTree.h"
+  #include "TObjArray.h"
+  #include "TStyle.h"
+  #include "TCanvas.h"
+  #include "TLine.h"
+  #include "TText.h"
+  #include "TParticle.h"
+  #include "TStopwatch.h"
+
+  #include "AliRun.h"
+  #include "AliPDG.h"
+  #include "AliV0vertex.h"
+#endif
+
+struct GoodVertex {
+  Int_t nlab,plab;
+  Int_t code;
+  Float_t px,py,pz;
+  Float_t x,y,z;
+};
+Int_t good_vertices(GoodVertex *gt, Int_t max);
+
+Int_t AliV0Comparison(Int_t code=310) { //Lambda=3122, LambdaBar=-3122
+   cerr<<"Doing comparison...\n";
+
+   const Double_t V0window=0.05, V0width=0.015; 
+   Double_t V0mass=0.5;
+   switch (code) {
+   case kK0Short:    V0mass=0.497672; break;
+   case kLambda0:    V0mass=1.115683; break;
+   case kLambda0Bar: V0mass=1.115683; break;
+   default: cerr<<"Invalid PDG code !\n"; return 1;
+   }
+
+   TStopwatch timer;
+
+   /*** Load reconstructed vertices ***/
+   TFile *vf=TFile::Open("AliV0vertices.root");
+   if (!vf->IsOpen()) {cerr<<"Can't open AliV0vertices.root !\n"; return 2;}
+   TObjArray varray(1000);
+   TTree *vTree=(TTree*)vf->Get("TreeV");
+   TBranch *branch=vTree->GetBranch("vertices");
+   Int_t nentr=(Int_t)vTree->GetEntries();
+   for (Int_t i=0; i<nentr; i++) {
+       AliV0vertex *iovertex=new AliV0vertex; branch->SetAddress(&iovertex);
+       vTree->GetEvent(i);
+       varray.AddLast(iovertex);
+   }
+
+   /*** Check if the file with the "good" vertices exists ***/
+   GoodVertex gv[1000];
+   Int_t ngood=0;
+   ifstream in("good_vertices");
+   if (in) {
+      cerr<<"Reading good vertices...\n";
+      while (in>>gv[ngood].nlab>>gv[ngood].plab>>gv[ngood].code>>
+                 gv[ngood].px>>gv[ngood].py>>gv[ngood].pz>>
+                 gv[ngood].x >>gv[ngood].y >>gv[ngood].z) {
+         ngood++;
+         cerr<<ngood<<'\r';
+         if (ngood==1000) {
+            cerr<<"Too many good vertices !\n";
+            break;
+         }
+      }
+      if (!in.eof()) cerr<<"Read error (good_vertices) !\n";
+   } else {
+     /*** generate a file with the "good" vertices ***/
+      cerr<<"Marking good vertices (this will take a while)...\n";
+      ngood=good_vertices(gv,1000);
+      ofstream out("good_vertices");
+      if (out) {
+         for (Int_t ngd=0; ngd<ngood; ngd++)            
+           out<<gv[ngd].nlab<<' '<<gv[ngd].plab<<' '<<gv[ngd].code<<' '<<
+                 gv[ngd].px<<' '<<gv[ngd].py<<' '<<gv[ngd].pz<<' '<<
+                 gv[ngd].x <<' '<<gv[ngd].y <<' '<<gv[ngd].z <<endl;
+      } else cerr<<"Can not open file (good_vertices) !\n";
+      out.close();
+   }
+
+   vf->Close();
+
+
+   TH1F *hp=new TH1F("hp","Angular Resolution",30,-30.,30.); //phi resolution 
+   hp->SetXTitle("(mrad)"); hp->SetFillColor(2);
+   TH1F *hl=new TH1F("hl","Lambda Resolution",30,-30,30);
+   hl->SetXTitle("(mrad)"); hl->SetFillColor(1); hl->SetFillStyle(3013); 
+   TH1F *hpt=new TH1F("hpt","Relative Pt Resolution",30,-10.,10.); 
+   hpt->SetXTitle("(%)"); hpt->SetFillColor(2); 
+
+   TH1F *hx=new TH1F("hx","Position Resolution (X,Y)",30,-3.,3.); //x res. 
+   hx->SetXTitle("(mm)"); hx->SetFillColor(6);
+   TH1F *hy=new TH1F("hy","Position Resolution (Y)",30,-3.,3.);   //y res
+   hy->SetXTitle("(mm)"); hy->SetFillColor(1); hy->SetFillStyle(3013);
+   TH1F *hz=new TH1F("hz","Position Resolution (Z)",30,-3.,3.);   //z res. 
+   hz->SetXTitle("(mm)"); hz->SetFillColor(6);
+
+
+   Double_t pmin=0.2, pmax=4.2; Int_t nchan=20;
+   TH1F *hgood=new TH1F("hgood","Good Vertices",nchan,pmin,pmax);    
+   TH1F *hfound=new TH1F("hfound","Found Vertices",nchan,pmin,pmax);
+   TH1F *hfake=new TH1F("hfake","Fake Vertices",nchan,pmin,pmax);
+   TH1F *hg=new TH1F("hg","Efficiency for Good Vertices",nchan,pmin,pmax);
+   hg->SetLineColor(4); hg->SetLineWidth(2);
+   TH1F *hf=new TH1F("hf","Probability of Fake Vertices",nchan,pmin,pmax);
+   hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
+
+   Double_t mmin=V0mass-V0window, mmax=V0mass+V0window;
+   TH1F *v0s =new TH1F("v0s","V0s Effective Mass",40, mmin, mmax);
+   v0s->SetXTitle("(GeV/c)"); v0s->SetFillColor(6);
+
+
+   Double_t pxg=0.,pyg=0.,ptg=0.;
+   Int_t nlab=-1, plab=-1;
+   Int_t i;
+   for (i=0; i<nentr; i++) {
+       AliV0vertex *vertex=(AliV0vertex*)varray.UncheckedAt(i);
+       nlab=TMath::Abs(vertex->GetNlabel());
+       plab=TMath::Abs(vertex->GetPlabel());
+
+       vertex->ChangeMassHypothesis(code);
+
+       Double_t mass=vertex->GetEffMass();
+       v0s->Fill(mass);
+
+       Int_t j;
+       for (j=0; j<ngood; j++) {
+          if (gv[j].code != vertex->GetPdgCode()) continue;
+          if (gv[j].nlab == nlab)
+          if (gv[j].plab == plab) break;
+       }
+
+       if (TMath::Abs(mass-V0mass)>V0width) continue;
+
+       Double_t px,py,pz; vertex->GetPxPyPz(px,py,pz);
+       Double_t pt=TMath::Sqrt(px*px+py*py);
+
+       if (j==ngood) {
+          hfake->Fill(pt);
+          cerr<<"Fake vertex: ("<<nlab<<","<<plab<<")\n";
+          continue;
+       }
+
+       pxg=gv[j].px; pyg=gv[j].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
+       Double_t phig=TMath::ATan2(pyg,pxg), phi=TMath::ATan2(py,px);
+       Double_t lamg=TMath::ATan2(gv[j].pz,ptg), lam=TMath::ATan2(pz,pt);
+       hp->Fill((phi - phig)*1000.);
+       hl->Fill((lam - lamg)*1000.);
+       hpt->Fill((1/pt - 1/ptg)/(1/ptg)*100.);
+
+       Double_t x,y,z; vertex->GetXYZ(x,y,z);
+       hx->Fill((x-gv[j].x)*10);
+       hy->Fill((y-gv[j].y)*10);
+       hz->Fill((z-gv[j].z)*10);
+
+       hfound->Fill(ptg);
+       gv[j].nlab=-1;
+
+   }
+   for (i=0; i<ngood; i++) {
+      if (gv[i].code != code) continue;
+      pxg=gv[i].px; pyg=gv[i].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
+      hgood->Fill(ptg);
+      nlab=gv[i].nlab; plab=gv[i].plab;
+      if (nlab < 0) continue;
+      cerr<<"Vertex ("<<nlab<<','<<plab<<") has not been found !\n";
+   }
+
+   varray.Delete();
+
+   Stat_t ng=hgood->GetEntries();
+   Stat_t nf=hfound->GetEntries();
+
+   cerr<<"Number of found vertices: "<<nentr<<" ("<<nf<<" in the peak)\n";
+   cerr<<"Number of good vertices: "<<ng<<endl;
+
+   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,580,610);
+   c1->Divide(2,2);
+
+   c1->cd(1);
+   gPad->SetFillColor(42); gPad->SetFrameFillColor(10); 
+   //hp->Fit("gaus");
+   hp->Draw();
+   hl->Draw("same"); c1->cd();
+
+   c1->cd(2);
+   gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
+   //hpt->Fit("gaus"); c1->cd();
+   hpt->Draw(); c1->cd();
+
+   c1->cd(3);
+   gPad->SetFillColor(42); gPad->SetFrameFillColor(10); 
+   //hx->Fit("gaus"); 
+   hx->Draw(); 
+   hy->Draw("same"); c1->cd();
+
+   c1->cd(4);
+   gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
+   //hz->Fit("gaus");
+   hz->Draw();
+
+
+   TCanvas *c2=new TCanvas("c2","",600,0,580,610);
+   c2->Divide(1,2);
+
+   c2->cd(1);
+   gPad->SetFillColor(42); gPad->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("V0 reconstruction efficiency");
+   hg->SetXTitle("Pt (GeV/c)");
+   hg->Draw();
+
+   TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
+   line1->Draw("same");
+   TLine *line2 = new TLine(pmin,0.9,pmax,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 vertices");
+   text->SetTextSize(0.05);
+   text->Draw();
+   text = new TText(0.453919,1.11408,"Good vertices");
+   text->SetTextSize(0.05);
+   text->Draw();
+
+
+   c2->cd(2);
+   gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
+   v0s->SetXTitle("(GeV/c)"); 
+   v0s->Fit("gaus","","",V0mass-V0width,V0mass+V0width);
+   Double_t max=v0s->GetMaximum();
+   TLine *line3 = new TLine(V0mass-V0width,0.,V0mass-V0width,max);
+   line3->Draw("same");
+   TLine *line4 = new TLine(V0mass+V0width,0.,V0mass+V0width,max);
+   line4->Draw("same");
+
+   timer.Stop(); timer.Print();
+
+   return 0;
+}
+
+Int_t good_vertices(GoodVertex *gv, Int_t max) {
+   Int_t nv=0;
+   /*** Get information about the cuts ***/
+   Double_t r2min=0.9*0.9;
+   Double_t r2max=2.9*2.9;
+
+   /*** Get labels of the "good" tracks ***/
+   Double_t dd; Int_t id, label[15000], ngt=0;
+   ifstream in("good_tracks_its");
+   if (!in) {
+     cerr<<"Can't open the file good_tracks_its \n";
+     return nv;
+   }
+   while (in>>label[ngt]>>id>>dd>>dd>>dd>>dd>>dd>>dd) {
+     ngt++;
+     if (ngt>=15000) {
+       cerr<<"Too many good ITS tracks !\n";
+       return nv;
+     }
+   }   
+   if (!in.eof()) {
+      cerr<<"Read error (good_tracks_its) !\n";
+      return nv;
+   }
+
+   /*** Get an access to the kinematics ***/
+   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 has not been found on galice.root !\n";
+     exit(5);
+   }
+
+   Int_t np=gAlice->GetEvent(0);
+   while (np--) {
+      cerr<<np<<'\r';
+      TParticle *p0=gAlice->Particle(np);
+
+      /*** only these V0s are "good" ***/
+      Int_t code=p0->GetPdgCode();
+      if (code!=kK0Short)
+      if (code!=kLambda0)
+      if (code!=kLambda0Bar) continue; 
+
+      /*** daughter tracks should be "good" ***/
+      Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
+      if (nlab==plab) continue;
+      if (nlab<0) continue;
+      if (plab<0) continue;
+      Int_t i;
+      for (i=0; i<ngt; i++) if (label[i]==nlab) break;
+      if (i==ngt) continue;
+      for (i=0; i<ngt; i++) if (label[i]==plab) break;
+      if (i==ngt) continue;
+
+      /*** fiducial volume ***/
+      TParticle *p=gAlice->Particle(nlab);
+      Double_t x=p->Vx(), y=p->Vy(), z=p->Vz(), r2=x*x+y*y;
+      if (r2<r2min) continue;
+      if (r2>r2max) continue;
+      gv[nv].code=code;
+      gv[nv].plab=plab; gv[nv].nlab=nlab;
+      gv[nv].px=p0->Px(); gv[nv].py=p0->Py(); gv[nv].pz=p0->Pz();
+      gv[nv].x=x; gv[nv].y=y; gv[nv].z=z;
+      nv++;
+   }
+   delete gAlice; gAlice=0;
+
+   file->Close();  
+
+   return nv;
+}
diff --git a/ITS/AliV0FindVertices.C b/ITS/AliV0FindVertices.C
new file mode 100644 (file)
index 0000000..fa10eb1
--- /dev/null
@@ -0,0 +1,35 @@
+#ifndef __CINT__
+  #include <iostream.h>
+  #include "AliV0vertexer.h"
+  #include "TFile.h"
+  #include "TStopwatch.h"
+#endif
+
+Int_t AliV0FindVertices() {
+   cerr<<"Looking for V0 vertices...\n";
+
+   TFile *out=TFile::Open("AliV0vertices.root","new");
+   if (!out->IsOpen()) {cerr<<"Delete old AliV0vertices.root !\n"; return 1;}
+
+   TFile *in=TFile::Open("AliITStracksV2.root");
+   if (!in->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 2;}
+
+   Double_t cuts[]={33.,  // max. allowed chi2
+                    0.015,// min. allowed negative daughter's impact parameter 
+                    0.015,// min. allowed positive daughter's impact parameter 
+                    0.060,// max. allowed DCA between the daughter tracks
+                    0.997,// max. allowed cosine of V0's pointing angle
+                    0.9,  // min. radius of the fiducial volume
+                    2.9   // max. radius of the fiducial volume
+                   };
+   TStopwatch timer;
+   AliV0vertexer *vertexer=new AliV0vertexer(cuts);
+   Int_t rc=vertexer->Tracks2V0vertices(in,out);
+   delete vertexer;
+   timer.Stop(); timer.Print();
+    
+   in->Close();
+   out->Close();
+
+   return rc;
+}
diff --git a/ITS/AliV0vertex.cxx b/ITS/AliV0vertex.cxx
new file mode 100644 (file)
index 0000000..ca305e3
--- /dev/null
@@ -0,0 +1,151 @@
+/**************************************************************************
+ * 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 V0 vertex class
+//
+//     Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
+//-------------------------------------------------------------------------
+#include <iostream.h>
+#include <TMath.h>
+
+#include "AliV0vertex.h"
+#include "AliITStrackV2.h"
+
+ClassImp(AliV0vertex)
+
+AliV0vertex::AliV0vertex() : TObject() {
+  //--------------------------------------------------------------------
+  // Default constructor  (K0s)
+  //--------------------------------------------------------------------
+  fPdgCode=kK0Short;
+  fEffMass=0.497672;
+  fChi2=1.e+33;
+  fPos[0]=fPos[1]=fPos[2]=0.;
+  fPosCov[0]=fPosCov[1]=fPosCov[2]=fPosCov[3]=fPosCov[4]=fPosCov[5]=0.;
+}
+
+AliV0vertex::AliV0vertex(const AliITStrackV2 &n, const AliITStrackV2 &p) {
+  //--------------------------------------------------------------------
+  // Main constructor
+  //--------------------------------------------------------------------
+  fPdgCode=kK0Short;
+  fNlab=n.GetLabel(); fPlab=p.GetLabel();
+
+  //Trivial estimation of the vertex parameters
+  Double_t pt, phi, x, par[5];
+  Double_t alpha, cs, sn;
+
+  n.GetExternalParameters(x,par); alpha=n.GetAlpha();
+  pt=1./TMath::Abs(par[4]);
+  phi=TMath::ASin(par[2]) + alpha;  
+  Double_t px1=pt*TMath::Cos(phi), py1=pt*TMath::Sin(phi), pz1=pt*par[3];
+  cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+  Double_t x1=x*cs - par[0]*sn;
+  Double_t y1=x*sn + par[0]*cs;
+  Double_t z1=par[1];
+  Double_t sx1=sn*sn*n.GetSigmaY2(), sy1=cs*cs*n.GetSigmaY2(); 
+
+  p.GetExternalParameters(x,par); alpha=p.GetAlpha();
+  pt=1./TMath::Abs(par[4]);
+  phi=TMath::ASin(par[2]) + alpha;  
+  Double_t px2=pt*TMath::Cos(phi), py2=pt*TMath::Sin(phi), pz2=pt*par[3];
+  cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+  Double_t x2=x*cs - par[0]*sn;
+  Double_t y2=x*sn + par[0]*cs;
+  Double_t z2=par[1];
+  Double_t sx2=sn*sn*p.GetSigmaY2(), sy2=cs*cs*p.GetSigmaY2(); 
+    
+  Double_t sz1=n.GetSigmaZ2(), sz2=p.GetSigmaZ2();
+  Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
+  Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
+  Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
+  fPos[0]=wx1*x1 + wx2*x2; fPos[1]=wy1*y1 + wy2*y2; fPos[2]=wz1*z1 + wz2*z2;
+
+  //fPos[0]=0.5*(x1+x2); fPos[1]=0.5*(y1+y2); fPos[2]=0.5*(z1+z2);
+  fNmom[0]=px1; fNmom[1]=py1; fNmom[2]=pz1; 
+  fPmom[0]=px2; fPmom[1]=py2; fPmom[2]=pz2;
+
+  Double_t e1=TMath::Sqrt(0.13957*0.13957 + px1*px1 + py1*py1 + pz1*pz1);
+  Double_t e2=TMath::Sqrt(0.13957*0.13957 + px2*px2 + py2*py2 + pz2*pz2);
+  fEffMass=TMath::Sqrt((e1+e2)*(e1+e2)-
+    (px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2));
+
+  fChi2=7.;   
+}
+
+void AliV0vertex::ChangeMassHypothesis(Int_t code) {
+  //--------------------------------------------------------------------
+  // This function changes the mass hypothesis for this V0
+  //--------------------------------------------------------------------
+  Double_t nmass, pmass;
+
+  switch (code) {
+  case kLambda0:
+    nmass=0.13957; pmass=0.93827; break;
+  case kLambda0Bar:
+    pmass=0.13957; nmass=0.93827; break;
+  case kK0Short: 
+    nmass=0.13957; pmass=0.13957; break;
+  default:
+    cerr<<"AliV0vertex::ChangeMassHypothesis: ";
+    cerr<<"invalide PDG code ! Assuming K0s...\n";
+    nmass=0.13957; pmass=0.13957; break;
+  }
+
+  Double_t px1=fNmom[0], py1=fNmom[1], pz1=fNmom[2]; 
+  Double_t px2=fPmom[0], py2=fPmom[1], pz2=fPmom[2];
+
+  Double_t e1=TMath::Sqrt(nmass*nmass + px1*px1 + py1*py1 + pz1*pz1);
+  Double_t e2=TMath::Sqrt(pmass*pmass + px2*px2 + py2*py2 + pz2*pz2);
+  fEffMass=TMath::Sqrt((e1+e2)*(e1+e2)-
+    (px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2));
+  
+  fPdgCode=code;
+}
+
+void AliV0vertex::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
+  //--------------------------------------------------------------------
+  // This function returns V0's momentum (global)
+  //--------------------------------------------------------------------
+  px=fNmom[0]+fPmom[0]; 
+  py=fNmom[1]+fPmom[1]; 
+  pz=fNmom[2]+fPmom[2]; 
+}
+
+void AliV0vertex::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
+  //--------------------------------------------------------------------
+  // This function returns V0's position (global)
+  //--------------------------------------------------------------------
+  x=fPos[0]; 
+  y=fPos[1]; 
+  z=fPos[2]; 
+}
+
+Double_t AliV0vertex::GetD(Double_t x0, Double_t y0, Double_t z0) const {
+  //--------------------------------------------------------------------
+  // This function returns V0's impact parameter
+  //--------------------------------------------------------------------
+  Double_t x=fPos[0],y=fPos[1],z=fPos[2];
+  Double_t px=fNmom[0]+fPmom[0];
+  Double_t py=fNmom[1]+fPmom[1];
+  Double_t pz=fNmom[2]+fPmom[2];
+
+  Double_t dx=(y0-y)*pz - (z0-z)*py; 
+  Double_t dy=(x0-x)*pz - (z0-z)*px;
+  Double_t dz=(x0-x)*py - (y0-y)*px;
+  Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
+  return d;
+}
diff --git a/ITS/AliV0vertex.h b/ITS/AliV0vertex.h
new file mode 100644 (file)
index 0000000..3de575a
--- /dev/null
@@ -0,0 +1,54 @@
+#ifndef ALIV0VERTEX_H
+#define ALIV0VERTEX_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                          V0 Vertex Class
+//
+//    Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch 
+//-------------------------------------------------------------------------
+
+#include <TObject.h>
+#include "AliPDG.h"
+
+class AliITStrackV2;
+
+class AliV0vertex : public TObject {
+public:
+  AliV0vertex();
+  AliV0vertex(const AliITStrackV2 &neg, const AliITStrackV2 &pos);
+
+  void ChangeMassHypothesis(Int_t code=kLambda0); 
+
+  Int_t GetPdgCode() const {return fPdgCode;}
+  Double_t GetEffMass() const {return fEffMass;}
+  Double_t GetChi2() const {return fChi2;}
+  void GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const;
+  void GetXYZ(Double_t &x, Double_t &y, Double_t &z) const;
+  Double_t GetD(Double_t x0=0.,Double_t y0=0.,Double_t z0=0.) const;
+  Int_t GetNlabel() const {return fNlab;}
+  Int_t GetPlabel() const {return fPlab;}
+
+private: 
+  Int_t fPdgCode;           // reconstructed V0's type (PDG code)
+  Double_t fEffMass;        // reconstructed V0's effective mass
+  Double_t fChi2;           // V0's chi2 value
+  Double_t fPos[3];         // V0's position (global)
+  Double_t fPosCov[6];      // covariance matrix of the vertex position
+
+  Int_t fNlab;              // label of the negative daughter
+  Double_t fNmom[3];        // momentum of the negative daughter (global)
+  Double_t fNmomCov[6];     // covariance matrix of the negative daughter mom.
+
+  Int_t fPlab;              // label of the positive daughter
+  Double_t fPmom[3];        // momentum of the positive daughter (global)
+  Double_t fPmomCov[6];     // covariance matrix of the positive daughter mom.
+
+  ClassDef(AliV0vertex,1)   // reconstructed V0 vertex
+};
+
+#endif
+
+
diff --git a/ITS/AliV0vertexer.cxx b/ITS/AliV0vertexer.cxx
new file mode 100644 (file)
index 0000000..077fb22
--- /dev/null
@@ -0,0 +1,283 @@
+/**************************************************************************
+ * 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 V0 vertexer class
+//
+//     Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
+//-------------------------------------------------------------------------
+#include <TFile.h>
+#include <TTree.h>
+#include <TObjArray.h>
+#include <iostream.h>
+
+#include "AliV0vertex.h"
+#include "AliV0vertexer.h"
+#include "AliITStrackV2.h"
+
+ClassImp(AliV0vertexer)
+
+Int_t AliV0vertexer::Tracks2V0vertices(const TFile *inp, TFile *out) {
+  //--------------------------------------------------------------------
+  //This function reconstructs V0 vertices
+  //--------------------------------------------------------------------
+   TFile *in=(TFile*)inp;
+   TDirectory *savedir=gDirectory; 
+
+   if (!in->IsOpen()) {
+     cerr<<"AliV0vertexer::Tracks2V0vertices(): ";
+     cerr<<"file with ITS tracks has not been open !\n";
+     return 1;
+   }
+
+   if (!out->IsOpen()) {
+     cerr<<"AliV0vertexer::Tracks2V0vertices(): ";
+     cerr<<"file for V0 vertices has not been open !\n";
+     return 2;
+   }
+
+   in->cd();
+
+   TTree *trkTree=(TTree*)in->Get("TreeT_ITS_0");
+   TBranch *branch=trkTree->GetBranch("tracks");
+   Int_t nentr=(Int_t)trkTree->GetEntries();
+
+   TObjArray negtrks(nentr/2);
+   TObjArray postrks(nentr/2);
+
+   Int_t nneg=0, npos=0, nvtx=0;
+
+   Int_t i;
+   for (i=0; i<nentr; i++) {
+       AliITStrackV2 *iotrack=new AliITStrackV2;
+       branch->SetAddress(&iotrack);
+       trkTree->GetEvent(i);
+       if (iotrack->Get1Pt() > 0.) {nneg++; negtrks.AddLast(iotrack);}
+       else {npos++; postrks.AddLast(iotrack);}
+   }   
+
+
+   out->cd();
+   TTree vtxTree("TreeV","Tree with V0 vertices");
+   AliV0vertex *ioVertex=0;
+   vtxTree.Branch("vertices","AliV0vertex",&ioVertex,32000,0);
+
+
+   for (i=0; i<nneg; i++) {
+      if (i%10==0) cerr<<nneg-i<<'\r';
+      AliITStrackV2 *ntrk=(AliITStrackV2 *)negtrks.UncheckedAt(i);
+      for (Int_t k=0; k<npos; k++) {
+         AliITStrackV2 *ptrk=(AliITStrackV2 *)postrks.UncheckedAt(k);
+
+         if (TMath::Abs(ntrk->GetD())<fDNmin) continue;
+         if (TMath::Abs(ptrk->GetD())<fDPmin) continue;
+
+         AliITStrackV2 nt(*ntrk), pt(*ptrk), *pnt=&nt, *ppt=&pt;
+
+         Double_t dca=PropagateToDCA(pnt,ppt);
+         if (dca > fDCAmax) continue;
+
+         AliV0vertex vertex(*pnt,*ppt);
+         if (vertex.GetChi2() > fChi2max) continue;
+        
+         /*  Think of something better here ! 
+         nt.PropagateToVertex(); if (TMath::Abs(nt.GetZ())<0.04) continue;
+         pt.PropagateToVertex(); if (TMath::Abs(pt.GetZ())<0.04) continue;
+        */
+
+         Double_t x,y,z; vertex.GetXYZ(x,y,z); 
+         Double_t r2=x*x + y*y; 
+         if (r2 > fRmax*fRmax) continue;
+         if (r2 < fRmin*fRmin) continue;
+
+         Double_t px,py,pz; vertex.GetPxPyPz(px,py,pz);
+         Double_t p2=px*px+py*py+pz*pz;
+         Double_t cost=(x*px+y*py+z*pz)/TMath::Sqrt(p2*(r2+z*z));
+         if (cost<fCPAmax) continue;
+
+         //vertex.ChangeMassHypothesis(); //default is Lambda0 
+
+         ioVertex=&vertex; vtxTree.Fill();
+
+         nvtx++;
+      }
+   }
+
+   cerr<<"Number of reconstructed V0 vertices: "<<nvtx<<endl;
+
+   vtxTree.Write();
+
+   negtrks.Delete();
+   postrks.Delete();
+
+   savedir->cd();
+   
+   delete trkTree;
+
+   return 0;
+}
+
+
+static void External2Helix(const AliITStrackV2 *t, Double_t helix[6]) { 
+  //--------------------------------------------------------------------
+  // External track parameters -> helix parameters 
+  //--------------------------------------------------------------------
+  Double_t alpha,x,cs,sn;
+  t->GetExternalParameters(x,helix); alpha=t->GetAlpha();
+
+  cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+  helix[5]=x*cs - helix[0]*sn;            // x0
+  helix[0]=x*sn + helix[0]*cs;            // y0
+//helix[1]=                               // z0
+  helix[2]=TMath::ASin(helix[2]) + alpha; // phi0
+//helix[3]=                               // tgl
+  helix[4]=helix[4]/t->GetConvConst();    // C
+}
+
+static void Evaluate(const Double_t *h, Double_t t,
+                     Double_t r[3],  //radius vector
+                     Double_t g[3],  //first defivatives
+                     Double_t gg[3]) //second derivatives
+{
+  //--------------------------------------------------------------------
+  // Calculate position of a point on a track and some derivatives
+  //--------------------------------------------------------------------
+  Double_t phase=h[4]*t+h[2];
+  Double_t sn=TMath::Sin(phase), cs=TMath::Cos(phase);
+
+  r[0] = h[5] + (sn - h[6])/h[4];
+  r[1] = h[0] - (cs - h[7])/h[4];  
+  r[2] = h[1] + h[3]*t;
+
+  g[0] = cs; g[1]=sn; g[2]=h[3];
+  
+  gg[0]=-h[4]*sn; gg[1]=h[4]*cs; gg[2]=0.;
+}
+
+Double_t AliV0vertexer::PropagateToDCA(AliITStrackV2 *n, AliITStrackV2 *p) {
+  //--------------------------------------------------------------------
+  // This function returns the DCA between two tracks
+  // The tracks will be moved to the point of DCA ! 
+  //--------------------------------------------------------------------
+  Double_t p1[8]; External2Helix(n,p1);
+  p1[6]=TMath::Sin(p1[2]); p1[7]=TMath::Cos(p1[2]);
+  Double_t p2[8]; External2Helix(p,p2);
+  p2[6]=TMath::Sin(p2[2]); p2[7]=TMath::Cos(p2[2]);
+
+
+  Double_t r1[3],g1[3],gg1[3]; Double_t t1=0.;
+  Evaluate(p1,t1,r1,g1,gg1);
+  Double_t r2[3],g2[3],gg2[3]; Double_t t2=0.;
+  Evaluate(p2,t2,r2,g2,gg2);
+
+  Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
+  Double_t dm=dx*dx + dy*dy + dz*dz;
+
+  Int_t max=27;
+  while (max--) {
+     Double_t gt1=-(dx*g1[0] + dy*g1[1] + dz*g1[2]);
+     Double_t gt2=+(dx*g2[0] + dy*g2[1] + dz*g2[2]);
+     Double_t h11=g1[0]*g1[0] - dx*gg1[0] + 
+                  g1[1]*g1[1] - dy*gg1[1] +
+                  g1[2]*g1[2] - dz*gg1[2];
+     Double_t h22=g2[0]*g2[0] + dx*gg2[0] + 
+                  g2[1]*g2[1] + dy*gg2[1] +
+                  g2[2]*g2[2] + dz*gg2[2];
+     Double_t h12=-(g1[0]*g2[0] + g1[1]*g2[1] + g1[2]*g2[2]);
+
+     Double_t det=h11*h22-h12*h12;
+
+     Double_t dt1,dt2;
+     if (TMath::Abs(det)<1.e-33) {
+        //(quasi)singular Hessian
+        dt1=-gt1; dt2=-gt2;
+     } else {
+        dt1=-(gt1*h22 - gt2*h12)/det; 
+        dt2=-(h11*gt2 - h12*gt1)/det;
+     }
+
+     if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
+
+     //check delta(phase1) ?
+     //check delta(phase2) ?
+
+     if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
+     if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
+        if ((gt1*gt1+gt2*gt2) > 1.e-4) 
+          cerr<<"AliV0vertexer::PropagateToDCA:"
+                 " stopped at not a stationary point !\n";
+        Double_t lmb=h11+h22; lmb=lmb-TMath::Sqrt(lmb*lmb-4*det);
+        if (lmb < 0.) 
+          cerr<<"AliV0vertexer::PropagateToDCA:"
+                 " stopped at not a minimum !\n";
+        break;
+     }
+
+     Double_t dd=dm;
+     for (Int_t div=1 ; ; div*=2) {
+        Evaluate(p1,t1+dt1,r1,g1,gg1);
+        Evaluate(p2,t2+dt2,r2,g2,gg2);
+        dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
+        dd=dx*dx + dy*dy + dz*dz;
+       if (dd<dm) break;
+        dt1*=0.5; dt2*=0.5;
+        if (div>512) {
+           cerr<<"AliV0vertexer::PropagateToDCA: overshoot !\n"; break;
+        }   
+     }
+     dm=dd;
+
+     t1+=dt1;
+     t2+=dt2;
+
+  }
+
+  if (max<=0) cerr<<"AliV0vertexer::PropagateToDCA: too many iterations !\n";  
+  
+  //propagate tracks to the points of DCA
+  Double_t cs=TMath::Cos(n->GetAlpha());
+  Double_t sn=TMath::Sin(n->GetAlpha());
+  Double_t x=r1[0]*cs + r1[1]*sn;
+  if (!n->PropagateTo(x,0.,0.)) {
+    //cerr<<"AliV0vertexer::PropagateToDCA: propagation failed !\n";
+    return 1.e+33;
+  }  
+
+  cs=TMath::Cos(p->GetAlpha());
+  sn=TMath::Sin(p->GetAlpha());
+  x=r2[0]*cs + r2[1]*sn;
+  if (!p->PropagateTo(x,0.,0.)) {
+    //cerr<<"AliV0vertexer::PropagateToDCA: propagation failed !\n";
+    return 1.e+33;
+  }  
+
+  return TMath::Sqrt(dm);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ITS/AliV0vertexer.h b/ITS/AliV0vertexer.h
new file mode 100644 (file)
index 0000000..014418c
--- /dev/null
@@ -0,0 +1,70 @@
+#ifndef ALIV0VERTEXER_H
+#define ALIV0VERTEXER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//------------------------------------------------------------------
+//                    V0 Vertexer Class
+//
+//   Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch 
+//------------------------------------------------------------------
+
+#include "TObject.h"
+
+class TFile;
+class AliITStrackV2;
+
+//_____________________________________________________________________________
+class AliV0vertexer : public TObject {
+public:
+  AliV0vertexer();
+  AliV0vertexer(const Double_t cuts[7]);
+  void SetCuts(const Double_t cuts[7]);
+
+  Int_t Tracks2V0vertices(const TFile *in, TFile *out);
+  Double_t PropagateToDCA(AliITStrackV2 *nt, AliITStrackV2 *pt);
+
+  void GetCuts(Double_t cuts[7]) const;
+
+private:
+  Double_t fChi2max;      // maximal allowed chi2 
+  Double_t fDNmin;        // min. allowed negative daughter's impact parameter
+  Double_t fDPmin;        // min. allowed positive daughter's impact parameter
+  Double_t fDCAmax;       // maximal allowed DCA between the daughter tracks 
+  Double_t fCPAmax;       // maximal allowed cosine of V0's pointing angle
+  Double_t fRmin, fRmax;  // max & min radii of the fiducial volume
+  
+  ClassDef(AliV0vertexer,1)  // V0 verterxer 
+};
+
+inline AliV0vertexer::AliV0vertexer() {
+ fChi2max=33.; 
+ fDNmin=0.015; fDPmin=0.015;
+ fDCAmax=0.01; fCPAmax=0.025; 
+ fRmin=0.5;    fRmax=2.5; 
+}
+
+inline AliV0vertexer::AliV0vertexer(const Double_t cuts[7]) {
+  fChi2max=cuts[0]; 
+  fDNmin=cuts[1];   fDPmin=cuts[2];
+  fDCAmax=cuts[3];  fCPAmax=cuts[4];
+  fRmin=cuts[5];    fRmax=cuts[6]; 
+}
+
+inline void AliV0vertexer::SetCuts(const Double_t cuts[7]) {
+  fChi2max=cuts[0]; 
+  fDNmin=cuts[1];   fDPmin=cuts[2];
+  fDCAmax=cuts[3];  fCPAmax=cuts[4];
+  fRmin=cuts[5];    fRmax=cuts[6]; 
+}
+
+inline void AliV0vertexer::GetCuts(Double_t cuts[7]) const {
+  cuts[0]=fChi2max; 
+  cuts[1]=fDNmin;   cuts[2]=fDPmin; 
+  cuts[3]=fDCAmax;  cuts[4]=fCPAmax;
+  cuts[5]=fRmin;    cuts[6]=fRmax; 
+}
+
+#endif
+
+
index 30af2ac..151fa16 100644 (file)
 #pragma link C++ class AliITSclusterV2+;
 #pragma link C++ class AliITStrackV2+;
 #pragma link C++ class AliITStrackerV2+;
-#pragma link C++ class  AliITSVertex+;
+#pragma link C++ class  AliV0vertex+;
+#pragma link C++ class  AliV0vertexer+;
 
+#pragma link C++ class  AliITSVertex+;
 #endif
index ce7913a..55bc2f0 100644 (file)
@@ -42,7 +42,8 @@ SRCS          = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \
                AliITSRad.cxx AliITSgeoinfo.cxx AliITSTrackerV1.cxx\
                AliITSvtest.cxx \
                 AliITSclusterV2.cxx AliITStrackV2.cxx AliITStrackerV2.cxx \
-               AliITSVertex.cxx                
+               AliITSVertex.cxx \
+               AliV0vertex.cxx AliV0vertexer.cxx
 #              AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \
 
 # Fortran sources
index 1b65068..d9ffec9 100644 (file)
@@ -31,7 +31,8 @@ SRCS          = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \
                AliITSRad.cxx AliITSgeoinfo.cxx AliITSTrackerV1.cxx\
                AliITSvtest.cxx \
                 AliITSclusterV2.cxx AliITStrackV2.cxx AliITStrackerV2.cxx \
-               AliITSVertex.cxx                
+               AliITSVertex.cxx \
+        AliV0vertex.cxx AliV0vertexer.cxx
 #              AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \
 
 HDRS:=  $(SRCS:.cxx=.h)
index 6fdcd8b..19dc66d 100644 (file)
@@ -16,14 +16,17 @@ class AliCluster;
 
 class AliKalmanTrack : public TObject {
 public:
-  AliKalmanTrack() { fLab=-3141593; fChi2=0; fN=0; }
-  AliKalmanTrack(const AliKalmanTrack &t) {fLab=t.fLab;fChi2=t.fChi2;fN=t.fN;}
+  AliKalmanTrack() { fLab=-3141593; fChi2=0; fN=0; fMass=0.13957;}
+  AliKalmanTrack(const AliKalmanTrack &t) {
+    fLab=t.fLab; fChi2=t.fChi2; fN=t.fN; fMass=t.fMass;
+  }
   virtual ~AliKalmanTrack(){};
   void SetLabel(Int_t lab) {fLab=lab;}
 
   Bool_t   IsSortable() const {return kTRUE;}
   Int_t    GetLabel()   const {return fLab;}
   Double_t GetChi2()    const {return fChi2;}
+  Double_t GetMass()    const {return fMass;}
   Int_t    GetNumberOfClusters() const {return fN;}
   virtual Int_t GetClusterIndex(Int_t i) const { //reserved for AliTracker
     printf("AliKalmanTrack::GetClusterIndex(Int_t i) must be overloaded !\n");
@@ -32,28 +35,26 @@ public:
 
   virtual Int_t Compare(const TObject *o) const {return 0;} 
 
-  virtual void GetExternalParameters(Double_t &xr, Double_t x[5]) const {}
-  virtual void GetExternalCovariance(Double_t cov[15]) const {}
+  virtual void GetExternalParameters(Double_t &xr, Double_t x[5]) const {;}
+  virtual void GetExternalCovariance(Double_t cov[15]) const {;}
 
-  virtual Double_t GetPredictedChi2(const AliCluster *cluster) const {return 0;}
+  virtual Double_t GetPredictedChi2(const AliCluster *cluster) const {return 0.;}
   virtual 
-  Int_t PropagateTo(Double_t xr,Double_t x0,Double_t rho,Double_t pm) {
-    return 0;
-  }
-  virtual Int_t Update(const AliCluster* c, Double_t chi2, UInt_t i) {
-    return 0;
-  }
+  Int_t PropagateTo(Double_t xr,Double_t x0,Double_t rho) {return 0;}
+  virtual Int_t Update(const AliCluster* c, Double_t chi2, UInt_t i) {return 0;}
 
   static void SetConvConst(Double_t cc) {fConvConst=cc;}
   Double_t GetConvConst() const {return fConvConst;}
 
 protected:
   void SetChi2(Double_t chi2) {fChi2=chi2;} 
+  void SetMass(Double_t mass) {fMass=mass;}
   void SetNumberOfClusters(Int_t n) {fN=n;} 
 
 private: 
   Int_t fLab;             // track label
   Double_t fChi2;         // total chi2 value for this track
+  Double_t fMass;         // mass hypothesis
   Int_t fN;               // number of associated clusters
 
   static Double_t fConvConst; //conversion constant cm -> GeV/c
index c2d7f01..8c3b66d 100644 (file)
 
 ClassImp(AliTracker)
 
+Double_t AliTracker::fX;
+Double_t AliTracker::fY;
+Double_t AliTracker::fZ;
+
 //__________________________________________________________________________
 void AliTracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
   //--------------------------------------------------------------------
index 0f13835..ce9769f 100644 (file)
@@ -16,15 +16,24 @@ class TFile;
 
 class AliTracker {
 public:
-  AliTracker(){}
+  AliTracker() { fX=fY=fZ=0.; }
   virtual ~AliTracker(){}
   virtual Int_t Clusters2Tracks(const TFile *in, TFile *out)=0;
   virtual Int_t PropagateBack(const TFile *in, TFile *out)=0;
+  static void SetVertex(Double_t *xyz) { fX=xyz[0]; fY=xyz[1]; fZ=xyz[2]; }
 
 //protected:
   virtual AliCluster *GetCluster(Int_t index) const=0;
   virtual void  UseClusters(const AliKalmanTrack *t, Int_t from=0) const;
   virtual void  CookLabel(AliKalmanTrack *t,Float_t wrong) const; 
+  Double_t GetX() const {return fX;}
+  Double_t GetY() const {return fY;}
+  Double_t GetZ() const {return fZ;}
+
+private:
+  static Double_t fX;  //X-coordinate of the primary vertex
+  static Double_t fY;  //Y-coordinate of the primary vertex
+  static Double_t fZ;  //Z-coordinate of the primary vertex
 
   ClassDef(AliTracker,1) //abstract tracker
 };
index 5bef125..056f67f 100644 (file)
@@ -1,87 +1,63 @@
+/****************************************************************************
+ *           Very important, delicate and rather obscure macro.             *
+ *                                                                          *
+ *               Creates list of "trackable" tracks,                        *
+ *             sorts tracks for matching with the ITS,                      *
+ *             calculates efficiency, resolutions etc.                      *
+ *                                                                          *
+ *           Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                 *
+ * with several nice improvements by: M.Ivanov, GSI, m.ivanov@gsi.de        *
+ ****************************************************************************/
+
 #ifndef __CINT__
   #include "alles.h"
   #include "AliTPCtracker.h"
 #endif
 
 struct GoodTrackTPC {
-  Int_t fEventN; //event number
   Int_t lab;
   Int_t code;
   Float_t px,py,pz;
   Float_t x,y,z;
 };
 
-Int_t AliTPCComparison(Int_t eventn=1) {
-   Int_t good_tracks_tpc(GoodTrackTPC *gt, Int_t max, Int_t eventn=1);
-
+Int_t AliTPCComparison(Int_t event=0) {
    cerr<<"Doing comparison...\n";
-   Int_t i;
-   gBenchmark->Start("AliTPCComparison");
-
-   TFile *cf=TFile::Open("AliTPCclusters.root");
-   if (!cf->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 1;}
-   AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
-   if (!digp) { cerr<<"TPC parameters have not been found !\n"; return 2; }
 
+   const Int_t MAX=20000;
+   Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const Int_t event);
 
+   gBenchmark->Start("AliTPCComparison");
 
-   ///////////
-   AliTPCtracker *tracker =0; 
-   TObjArray tarray(2000);
-   AliTPCtrack *iotrack=0;
-   Int_t nentr= 0;
-   Int_t eventptr[1000];
+   Int_t nentr=0,i=0; TObjArray tarray(MAX);
+   { /*Load tracks*/
    TFile *tf=TFile::Open("AliTPCtracks.root");
-   TTree *tracktree=0;
-   // Load clusters
-   eventptr[0]=0;
-   eventptr[1]=0;
-   char   tname[100];  //Y.B.
-   for (Int_t event=0;event<eventn; event++){
-     cf->cd();
-     tracker = new AliTPCtracker(digp,event);
-     tracker->LoadInnerSectors();
-     tracker->LoadOuterSectors();
-     //Y.B. char   tname[100];
-     if (eventn==-1) {
-       sprintf(tname,"TreeT_TPC");
-     }
-     else {
-       sprintf(tname,"TreeT_TPC_%d",event);
-     }
+   if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
+
+   char tname[100]; sprintf(tname,"TreeT_TPC_%d",event);
+   TTree *tracktree=(TTree*)tf->Get(tname);
+   if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n"; return 4;}
 
-     // Load tracks
-     if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
-     tracktree=(TTree*)tf->Get(tname);
-     if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n"; return 4;}
-     TBranch *tbranch=tracktree->GetBranch("tracks");
-     Int_t nentr0=(Int_t)tracktree->GetEntries();
-     nentr+=nentr0;
-     for (i=0; i<nentr0; i++) {
+   TBranch *tbranch=tracktree->GetBranch("tracks");
+   nentr=(Int_t)tracktree->GetEntries();
+   AliTPCtrack *iotrack=0;
+   for (i=0; i<nentr; i++) {
        iotrack=new AliTPCtrack;
        tbranch->SetAddress(&iotrack);
        tracktree->GetEvent(i);
-       tracker->CookLabel(iotrack,0.1);
        tarray.AddLast(iotrack);
-     }   
-     eventptr[event+1] = nentr;  //store end of the event
-     delete tracker;
-     delete tracktree; //Thanks to Mariana Bondila
-   }
+   }   
+   delete tracktree; //Thanks to Mariana Bondila
    tf->Close();
-   delete digp; //Thanks to Mariana Bondila
-   cf->Close();
-  
+   }
 
-/////////////////////////////////////////////////////////////////////////
-   const Int_t MAX=45000;
+   /* Generate a list of "good" tracks */
    GoodTrackTPC gt[MAX];
-
    Int_t ngood=0;
    ifstream in("good_tracks_tpc");
    if (in) {
       cerr<<"Reading good tracks...\n";
-      while (in>>gt[ngood].fEventN>>gt[ngood].lab>>gt[ngood].code>>
+      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++;
@@ -94,29 +70,18 @@ Int_t AliTPCComparison(Int_t eventn=1) {
       if (!in.eof()) cerr<<"Read error (good_tracks_tpc) !\n";
    } else {
       cerr<<"Marking good tracks (this will take a while)...\n";
-      ngood=good_tracks_tpc(gt,MAX,eventn);   //mi change
+      ngood=good_tracks_tpc(gt,MAX, event);
       ofstream out("good_tracks_tpc");
       if (out) {
          for (Int_t ngd=0; ngd<ngood; ngd++)            
-          out<<gt[ngd].fEventN<<' '<<gt[ngd].lab<<' '<<gt[ngd].code<<' '<<
+           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_tpc) !\n";
       out.close();
-
-      cerr<<"Preparing tracks for matching with the ITS...\n";
-      tarray.Sort();
-      tf=TFile::Open("AliTPCtracks.root","recreate");
-      tracktree=new TTree(tname,"Tree with \"forward\" TPC tracks");
-      tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
-      for (i=0; i<nentr; i++) {
-          iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
-          tracktree->Fill();
-      }
-      tracktree->Write();
-      tf->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);
@@ -125,54 +90,41 @@ Int_t AliTPCComparison(Int_t eventn=1) {
    TH1F *hmpt=new TH1F("hmpt","Relative Pt resolution (pt>4GeV/c)",30,-60,60); 
    hmpt->SetFillColor(6);
 
-   AliTPCtrack *trk=(AliTPCtrack*)tarray.UncheckedAt(0);
-   Double_t pmin=0.1*(100/0.299792458/0.2/trk->GetConvConst());
-   Double_t pmax=6.0+pmin;
-
-   TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);    
-   TH1F *hfound=new TH1F("hfound","Found tracks",30,pmin,pmax);
-   TH1F *hfake=new TH1F("hfake","Fake tracks",30,pmin,pmax);
-   TH1F *hg=new TH1F("hg","",30,pmin,pmax); //efficiency for good tracks
+   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,pmin,pmax);
+   TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,0.1,6.1);
    hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
 
    TH1F *he =new TH1F("he","dE/dX for pions with 0.4<p<0.5 GeV/c",50,0.,100.);
-   TH2F *hep=new TH2F("hep","dE/dX vs momentum",50,0.,2.,50,0.,250.);
+   TH2F *hep=new TH2F("hep","dE/dX vs momentum",50,0.,2.,50,0.,400.);
    hep->SetMarkerStyle(8);
    hep->SetMarkerSize(0.4);
 
-   Int_t mingood = ngood;  //MI change
-   Int_t * track_notfound = new Int_t[ngood];
-   Int_t   itrack_notfound =0;
-   Int_t * track_fake  = new Int_t[ngood];
-   Int_t   itrack_fake = 0;
-   Int_t * track_multifound = new Int_t[ngood];
-   Int_t * track_multifound_n = new Int_t[ngood];
-   Int_t   itrack_multifound =0;
+   //MI change
+   Int_t mingood=ngood;
+   Int_t track_notfound[MAX], itrack_notfound=0;
+   Int_t track_fake[MAX], itrack_fake=0;
+   Int_t track_multifound[MAX], track_multifound_n[MAX], itrack_multifound=0;
 
    while (ngood--) {
       Int_t lab=gt[ngood].lab,tlab=-1;
       Float_t ptg=
       TMath::Sqrt(gt[ngood].px*gt[ngood].px + gt[ngood].py*gt[ngood].py);
 
-      if (ptg<pmin) continue;
+      if (ptg<1e-33) continue; // for those not crossing 0 pad row
 
       hgood->Fill(ptg);
-      Int_t ievent = gt[ngood].fEventN;
 
       AliTPCtrack *track=0;
-      //      for (i=0; i<nentr; i++) {
-      for (i=eventptr[ievent]; i<eventptr[ievent+1]; i++) {
-
+      for (i=0; i<nentr; i++) {
           track=(AliTPCtrack*)tarray.UncheckedAt(i);
           tlab=track->GetLabel();
           if (lab==TMath::Abs(tlab)) break;
       }
-      //      if (i==nentr) {
-      if (i==eventptr[ievent+1]) {
-
-       //cerr<<"Track "<<lab<<" was not found !\n";
+      if (i==nentr) {
        track_notfound[itrack_notfound++]=lab;
         continue;
       }
@@ -181,29 +133,23 @@ Int_t AliTPCComparison(Int_t eventn=1) {
       Int_t micount=0;
       Int_t mi;
       AliTPCtrack * mitrack;
-      //      for (mi=0; mi<nentr; mi++) {
-      for (mi=eventptr[ievent]; mi<eventptr[ievent+1]; mi++) {
-
+      for (mi=0; mi<nentr; mi++) {
        mitrack=(AliTPCtrack*)tarray.UncheckedAt(mi);          
        if (lab==TMath::Abs(mitrack->GetLabel())) micount++;
       }
       if (micount>1) {
-       //cout<<"Track no. "<<lab<<" found "<<micount<<"  times\n";
        track_multifound[itrack_multifound]=lab;
        track_multifound_n[itrack_multifound]=micount;
        itrack_multifound++;
       }
       
-
-      //
-      Double_t xk=gt[ngood].x;//digp->GetPadRowRadii(0,0);
+      Double_t xk=gt[ngood].x;
       track->PropagateTo(xk);
 
       if (lab==tlab) hfound->Fill(ptg);
       else { 
        track_fake[itrack_fake++]=lab;
        hfake->Fill(ptg); 
-       //cerr<<lab<<" fake\n";
       }
 
       Double_t par[5]; track->GetExternalParameters(xk,par);
@@ -235,19 +181,7 @@ Int_t AliTPCComparison(Int_t eventn=1) {
 
    }
 
-   
-   Stat_t ng=hgood->GetEntries(); cerr<<"Good tracks "<<ng<<endl;
-   Stat_t nf=hfound->GetEntries();
-   if (ng!=0) 
-      cerr<<"\n\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
-   
-   //MI change  - addition
-   cout<<"Total number of found tracks ="<<nentr<<"\n";
-   cout<<"Total number of \"good\" tracks ="<<mingood<<"\n";
-
-
    cout<<"\nList of Not found tracks :\n";
-   //   Int_t i;
    for ( i = 0; i< itrack_notfound; i++){
      cout<<track_notfound[i]<<"\t";
      if ((i%5)==4) cout<<"\n";
@@ -259,11 +193,17 @@ Int_t AliTPCComparison(Int_t eventn=1) {
    }
     cout<<"\nList of multiple found tracks :\n";
    for ( i=0; i<itrack_multifound; i++) {
-       cout<<"id.   "<<track_multifound[i]<<"     found - "<<track_multifound_n[i]<<"times\n";
+       cout<<"id.   "<<track_multifound[i]
+           <<"     found - "<<track_multifound_n[i]<<"times\n";
    }
-   cout<<"\n\n\n";
 
-   
+   Stat_t ng=hgood->GetEntries(), nf=hfound->GetEntries();
+   if (ng!=0) cerr<<"\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
+   cout<<"Total number of found tracks ="<<nentr<<endl;
+   cout<<"Total number of \"good\" tracks ="
+       <<mingood<<"   (selected for comparison: "<<ng<<')'<<endl<<endl;
+
+
    gStyle->SetOptStat(111110);
    gStyle->SetOptFit(1);
 
@@ -295,9 +235,9 @@ Int_t AliTPCComparison(Int_t eventn=1) {
    hg->SetXTitle("Pt (GeV/c)");
    hg->Draw();
 
-   TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
+   TLine *line1 = new TLine(0.1,1.0,6.1,1.0); line1->SetLineStyle(4);
    line1->Draw("same");
-   TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
+   TLine *line2 = new TLine(0.1,0.9,6.1,0.9); line2->SetLineStyle(4);
    line2->Draw("same");
 
    hf->SetFillColor(1);
@@ -335,8 +275,8 @@ Int_t AliTPCComparison(Int_t eventn=1) {
 }
 
 
-Int_t good_tracks_tpc(GoodTrackTPC *gt, Int_t max, Int_t eventn) {
-  //eventn  - number of events in file
+Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const Int_t event) {
+   Int_t nt=0;
 
    TFile *file=TFile::Open("galice.root");
    if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
@@ -346,7 +286,7 @@ Int_t good_tracks_tpc(GoodTrackTPC *gt, Int_t max, Int_t eventn) {
      exit(5);
    }
 
-   //   Int_t np=gAlice->GetEvent(0); //MI change   
+   Int_t np=gAlice->GetEvent(event);   
 
    AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
    Int_t ver = TPC->IsVersion(); 
@@ -359,151 +299,171 @@ Int_t good_tracks_tpc(GoodTrackTPC *gt, Int_t max, Int_t eventn) {
    Int_t nrow_up=digp->GetNRowUp();
    Int_t nrows=digp->GetNRowLow()+nrow_up;
    Int_t zero=digp->GetZeroSup();
-   Int_t gap=Int_t(0.125*nrows);
+   Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
    Int_t good_number=Int_t(0.4*nrows);
 
-
-   //MI change
-   Int_t nt=0;  //reset counter
-   char treeName[100];  //declare event identifier
-       
-   for (Int_t event=0;event<eventn;event++){
-
-     Int_t np=gAlice->GetEvent(event);   
-     Int_t *good=new Int_t[np];
-     for (Int_t ii=0; ii<np; ii++) good[ii]=0;
-     
-     
-     //MI change to be possible compile macro
-     //definition out of the swith statemnet
-     Int_t sectors_by_rows=0;
-     TTree *TD=0;
-     AliSimDigits da, *digits=&da;
-     Int_t *count=0;
-     switch (ver) {
-     case 1:
-       {
-        TFile *cf=TFile::Open("AliTPCclusters.root");
-        if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
-        AliTPCClustersArray *ca=new AliTPCClustersArray;
-        ca->Setup(digp);
-        ca->SetClusterType("AliTPCcluster");
-        //ca->ConnectTree("Segment Tree");
-        ca->ConnectTree("TreeC_TPC_0");
-        Int_t nrows=Int_t(ca->GetTree()->GetEntries());
-        for (Int_t n=0; n<nrows; n++) {
-          AliSegmentID *s=ca->LoadEntry(n);
-          Int_t sec,row;
-          digp->AdjustSectorRow(s->GetID(),sec,row);
-          AliTPCClustersRow &clrow = *ca->GetRow(sec,row);
-          Int_t ncl=clrow.GetArray()->GetEntriesFast();
-          while (ncl--) {
-            AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
-            Int_t lab=c->GetLabel(0);
-            if (lab<0) continue; //noise cluster
-            lab=TMath::Abs(lab);
-            if (sec>=digp->GetNInnerSector())
-              if (row==nrow_up-1    ) good[lab]|=0x1000;
-            if (sec>=digp->GetNInnerSector())
-              if (row==nrow_up-1-gap) good[lab]|=0x800;
-            good[lab]++;
-          }
-          ca->ClearRow(sec,row);
-        }
-        cf->Close();
-       }
+   Int_t *good=new Int_t[np];
+   Int_t i;
+   for (i=0; i<np; i++) good[i]=0;
+
+
+   switch (ver) {
+   case 1:
+     {
+      char cname[100]; sprintf(cname,"TreeC_TPC_%d",event);
+      TFile *cf=TFile::Open("AliTPCclusters.root");
+      if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
+      AliTPCClustersArray *pca=new AliTPCClustersArray, &ca=*pca;
+      ca.Setup(digp);
+      ca.SetClusterType("AliTPCcluster");
+      ca.ConnectTree(cname);
+      Int_t nrows=Int_t(ca.GetTree()->GetEntries());
+      for (Int_t n=0; n<nrows; n++) {
+          AliSegmentID *s=ca.LoadEntry(n);
+          Int_t sec,row;
+          digp->AdjustSectorRow(s->GetID(),sec,row);
+          AliTPCClustersRow &clrow = *ca.GetRow(sec,row);
+          Int_t ncl=clrow.GetArray()->GetEntriesFast();
+          while (ncl--) {
+              AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
+              Int_t lab=c->GetLabel(0);
+              if (lab<0) continue; //noise cluster
+              lab=TMath::Abs(lab);
+
+              if (sec>=digp->GetNInnerSector())
+                 if (row==nrow_up-1) good[lab]|=0x4000;
+              if (sec>=digp->GetNInnerSector())
+                 if (row==nrow_up-1-gap) good[lab]|=0x1000;
+
+              if (sec>=digp->GetNInnerSector())
+                 if (row==nrow_up-1-shift) good[lab]|=0x2000;
+              if (sec>=digp->GetNInnerSector())
+                 if (row==nrow_up-1-gap-shift) good[lab]|=0x800;
+
+              good[lab]++;
+          }
+          ca.ClearRow(sec,row);
+      }
+      delete pca;
+      cf->Close();
+     }
      break;
-     case 2:
-
-       sprintf(treeName,"TreeD_75x40_100x60_%d",event);       
-       TD=(TTree*)gDirectory->Get(treeName);
-       TD->GetBranch("Segment")->SetAddress(&digits);
-       count = new Int_t[np];
-       Int_t i;
-       for (i=0; i<np; i++) count[i]=0;
-       sectors_by_rows=(Int_t)TD->GetEntries();
-       for (i=0; i<sectors_by_rows; i++) {
-        if (!TD->GetEvent(i)) continue;
-        Int_t sec,row;
-        digp->AdjustSectorRow(digits->GetID(),sec,row);
-        cerr<<sec<<' '<<row<<"                                     \r";
-        digits->First();
-        do { //Many thanks to J.Chudoba who noticed this
-          Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
-          Short_t dig = digits->GetDigit(it,ip);
-          Int_t idx0=digits->GetTrackID(it,ip,0); 
-          Int_t idx1=digits->GetTrackID(it,ip,1);
-          Int_t idx2=digits->GetTrackID(it,ip,2);
-          if (idx0>=0 && dig>=zero) count[idx0]+=1;
-          if (idx1>=0 && dig>=zero) count[idx1]+=1;
-          if (idx2>=0 && dig>=zero) count[idx2]+=1;
+   case 2:
+     {
+      char dname[100]; sprintf(dname,"TreeD_75x40_100x60_%d",event);
+      TTree *TD=(TTree*)gDirectory->Get(dname);
+      AliSimDigits da, *digits=&da;
+      TD->GetBranch("Segment")->SetAddress(&digits);
+
+      Int_t *count = new Int_t[np];
+      Int_t i;
+      for (i=0; i<np; i++) count[i]=0;
+
+      Int_t sectors_by_rows=(Int_t)TD->GetEntries();
+      for (i=0; i<sectors_by_rows; i++) {
+          if (!TD->GetEvent(i)) continue;
+          Int_t sec,row;
+          digp->AdjustSectorRow(digits->GetID(),sec,row);
+          cerr<<sec<<' '<<row<<"                                     \r";
+          digits->First();
+          do { //Many thanks to J.Chudoba who noticed this
+              Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
+              Short_t dig = digits->GetDigit(it,ip);
+              Int_t idx0=digits->GetTrackID(it,ip,0); 
+              Int_t idx1=digits->GetTrackID(it,ip,1);
+              Int_t idx2=digits->GetTrackID(it,ip,2);
+              if (idx0>=0 && dig>=zero) count[idx0]+=1;
+              if (idx1>=0 && dig>=zero) count[idx1]+=1;
+              if (idx2>=0 && dig>=zero) count[idx2]+=1;
           } while (digits->Next());
-        for (Int_t j=0; j<np; j++) {
-          if (count[j]>1) {
-            if (sec>=digp->GetNInnerSector())
-              if (row==nrow_up-1    ) good[j]|=0x1000;
-            if (sec>=digp->GetNInnerSector())
-              if (row==nrow_up-1-gap) good[j]|=0x800;
-            good[j]++;
-          }
-          count[j]=0;
-        }
-       }
-       delete[] count;
-       delete TD; //Thanks to Mariana Bondila
-       break;
-     default:
-       cerr<<"Invalid TPC version !\n";
-       file->Close();
-       exit(7);
-     }
-     
-     TTree *TH=gAlice->TreeH();
-     Int_t npart=(Int_t)TH->GetEntries();
-
-     while (npart--) {
-       AliTPChit *hit0=0;
-       
-       TPC->ResetHits();
-       TH->GetEvent(npart);
-       AliTPChit * hit = (AliTPChit*) TPC->FirstHit(-1);
-       while (hit){
-        if (hit->fQ==0.) break;
-       hit =  (AliTPChit*) TPC->NextHit();
-       }
-       if (hit) {
-        hit0 = new AliTPChit(*hit); //Make copy of hit
-        hit = hit0;
-       }
-       else continue;
-       AliTPChit *hit1=(AliTPChit*)TPC->NextHit();     
-       if (hit1==0) continue;
-       if (hit1->fQ != 0.) continue;
-       Int_t i=hit->Track();
-       TParticle *p = (TParticle*)gAlice->Particle(i);
-       
-       if (p->GetFirstMother()>=0) continue;  //secondary particle
-       if (good[i] < 0x1000+0x800+2+good_number) continue;
-       if (p->Pt()<0.100) continue;
-       if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
-       
-       gt[nt].lab=i;
-       gt[nt].code=p->GetPdgCode();
-       //**** px py pz - in global coordinate system, x y z - in local !
-       gt[nt].px=hit->X(); gt[nt].py=hit->Y(); gt[nt].pz=hit->Z();
-       Float_t cs,sn; digp->AdjustCosSin(hit1->fSector,cs,sn);
-       gt[nt].fEventN=event;  //MI change
-       gt[nt].x = hit1->X()*cs + hit1->Y()*sn;
-       gt[nt].y =-hit1->X()*sn + hit1->Y()*cs;
-       gt[nt].z = hit1->Z();
-       nt++;   
-       if (hit0) delete hit0;
-       cerr<<i<<"                \r";
-       if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
+          for (Int_t j=0; j<np; j++) {
+              if (count[j]>1) {
+                 if (sec>=digp->GetNInnerSector())
+                  if (row==nrow_up-1    ) good[j]|=0x4000;
+                 if (sec>=digp->GetNInnerSector())
+                  if (row==nrow_up-1-gap) good[j]|=0x1000;
+
+                 if (sec>=digp->GetNInnerSector())
+                  if (row==nrow_up-1-shift) good[j]|=0x2000;
+                 if (sec>=digp->GetNInnerSector())
+                  if (row==nrow_up-1-gap-shift) good[j]|=0x800;
+                 good[j]++;
+              }
+              count[j]=0;
+          }
+      }
+      delete[] count;
+      delete TD; //Thanks to Mariana Bondila
      }
-     delete[] good;
+      break;
+   default:
+      cerr<<"Invalid TPC version !\n";
+      file->Close();
+      exit(7);
    }
+
+
+   /** select tracks which are "good" enough **/
+   for (i=0; i<np; i++) {
+      if ((good[i]&0x5000) != 0x5000)
+      if ((good[i]&0x2800) != 0x2800) continue;
+      if ((good[i]&0x7FF ) < good_number) continue;
+
+      TParticle *p = (TParticle*)gAlice->Particle(i);
+
+      if (p->Pt()<0.100) continue;
+      if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
+
+      Int_t j=p->GetFirstMother();
+      if (j>=0) {
+         TParticle *pp = (TParticle*)gAlice->Particle(j);
+         if (pp->GetFirstMother()>=0) continue;//only one decay is allowed
+      }
+
+      gt[nt].lab=i;
+      gt[nt].code=p->GetPdgCode();
+      gt[nt].px=0.; gt[nt].py=0.; gt[nt].pz=0.;
+      gt[nt].x=0.; gt[nt].z=0.; gt[nt].z=0.;
+      nt++;
+      if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
+      cerr<<np-i<<"                \r";
+   }
+
+   /** check if there is also information at the entrance of the TPC **/
+   TTree *TH=gAlice->TreeH(); np=(Int_t)TH->GetEntries();
+   for (i=0; i<np; i++) {
+      TPC->ResetHits();
+      TH->GetEvent(i);
+      AliTPChit *phit = (AliTPChit*)TPC->FirstHit(-1);
+      for ( ; phit; phit=(AliTPChit*)TPC->NextHit() ) {
+        if (phit->fQ !=0. ) continue;
+
+         Double_t px=phit->X(), py=phit->Y(), pz=phit->Z();
+
+         if ((phit=(AliTPChit*)TPC->NextHit())==0) break;
+         if (phit->fQ != 0.) continue;
+
+         Double_t x=phit->X(), y=phit->Y(), z=phit->Z();
+        if (TMath::Sqrt(x*x+y*y)>90.) continue;
+
+         Int_t j, lab=phit->Track();
+         for (j=0; j<nt; j++) {if (gt[j].lab==lab) break;}
+         if (j==nt) continue;         
+
+        // (px,py,pz) - in global coordinate system, (x,y,z) - in local !
+         gt[j].px=px; gt[j].py=py; gt[j].pz=pz;
+         Float_t cs,sn; digp->AdjustCosSin(phit->fSector,cs,sn);
+         gt[j].x = x*cs + y*sn;
+         gt[j].y =-x*sn + y*cs;
+         gt[j].z = z;
+      }
+      cerr<<np-i<<"                \r";
+   }
+
+   delete[] good;
+
    delete gAlice; gAlice=0;
+
    file->Close();
    gBenchmark->Stop("AliTPCComparison");
    gBenchmark->Show("AliTPCComparison");   
index 7948dd3..a9ea84d 100644 (file)
@@ -1,5 +1,6 @@
 #ifndef __CINT__
   #include <iostream.h>
+  #include "AliTPCParam.h"
   #include "AliTPCtracker.h"
 
   #include "TFile.h"
@@ -20,10 +21,12 @@ Int_t AliTPCFindTracks(Int_t eventn=1) {
  
    TStopwatch timer;
 
+   Int_t rc=0;
    for (Int_t i=0;i<eventn;i++){
      printf("Processing event %d\n",i);
      AliTPCtracker *tracker = new AliTPCtracker(par,i);
-     Int_t rc=tracker->Clusters2Tracks(0,out);
+     //Double_t xyz[]={0.,0.,0.}; tracker->SetVertex(xyz); //primary vertex
+     rc=tracker->Clusters2Tracks(0,out);
      delete tracker;
    }
    timer.Stop(); timer.Print();
@@ -33,5 +36,5 @@ Int_t AliTPCFindTracks(Int_t eventn=1) {
    in->Close();
    out->Close();
 
-   return 0;
+   return rc;
 }
index 36fd9f6..b2fc940 100644 (file)
@@ -60,6 +60,7 @@ AliTPCtrack::AliTPCtrack(const AliKalmanTrack& t,Double_t alpha) {
   //-----------------------------------------------------------------
   SetLabel(t.GetLabel());
   SetChi2(0.);
+  SetMass(t.GetMass());
   SetNumberOfClusters(0);
   //SetConvConst(t.GetConvConst());
 
@@ -170,8 +171,7 @@ Double_t AliTPCtrack::GetPredictedChi2(const AliCluster *c) const
 }
 
 //_____________________________________________________________________________
-Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
-{
+Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho) {
   //-----------------------------------------------------------------
   // This function propagates a track to a reference plane x=xk.
   //-----------------------------------------------------------------
@@ -223,7 +223,7 @@ Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
   //Multiple scattering******************
   Double_t d=sqrt((x1-fX)*(x1-fX)+(y1-fP0)*(y1-fP0)+(z1-fP1)*(z1-fP1));
   Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
-  Double_t beta2=p2/(p2 + pm*pm);
+  Double_t beta2=p2/(p2 + GetMass()*GetMass());
   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
   //Double_t theta2=1.0259e-6*10*10/20/(beta2*p2)*d*rho;
 
@@ -241,14 +241,14 @@ Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
   Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
   if (x1 < x2) dE=-dE;
   cc=fP4;
-  fP4*=(1.- sqrt(p2+pm*pm)/p2*dE);
+  fP4*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
   fP2+=fX*(fP4-cc);
 
   return 1;
 }
 
 //_____________________________________________________________________________
-Int_t AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm) 
+Int_t AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho) 
 {
   //-----------------------------------------------------------------
   // This function propagates tracks to the "vertex".
@@ -257,7 +257,7 @@ Int_t AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm)
   Double_t tgf=-fP2/(fP4*fP0 + sqrt(1-c*c));
   Double_t snf=tgf/sqrt(1.+ tgf*tgf);
   Double_t xv=(fP2+snf)/fP4;
-  return PropagateTo(xv,x0,rho,pm);
+  return PropagateTo(xv,x0,rho);
 }
 
 //_____________________________________________________________________________
index 1626f60..8139b09 100644 (file)
@@ -39,10 +39,9 @@ public:
               const Double_t cc[15], Double_t xr, Double_t alpha); 
   AliTPCtrack(const AliKalmanTrack& t, Double_t alpha);
   AliTPCtrack(const AliTPCtrack& t);
-  Int_t PropagateToVertex(
-                    Double_t x0=36.66,Double_t rho=1.2e-3,Double_t pm=0.139);
+  Int_t PropagateToVertex(Double_t x0=36.66,Double_t rho=1.2e-3);
   Int_t Rotate(Double_t angle);
-  void SetdEdx(Float_t dedx) {fdEdx=dedx;}
+  void SetdEdx(Double_t dedx) {fdEdx=dedx;}
 
   Double_t GetX()     const {return fX;}
   Double_t GetAlpha() const {return fAlpha;}
@@ -77,8 +76,7 @@ public:
 //****************************************************** 
 
   virtual Double_t GetPredictedChi2(const AliCluster *cluster) const;
-  Int_t PropagateTo(Double_t xr,
-                    Double_t x0=28.94,Double_t rho=0.9e-3,Double_t pm=0.139);
+  Int_t PropagateTo(Double_t xr,Double_t x0=28.94,Double_t rho=0.9e-3);
   Int_t Update(const AliCluster* c, Double_t chi2, UInt_t i);
   void ResetCovariance();
 
@@ -87,7 +85,7 @@ private:
   Double_t fAlpha;          // Rotation angle the local (TPC sector)
                             // coordinate system and the global ALICE one.
 
-  Double_t fdEdx;           // dE/dx 
+  Double_t fdEdx;           // dE/dx
 
   Double_t fP0;             // Y-coordinate of a track
   Double_t fP1;             // Z-coordinate of a track
index 7538292..c706776 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.14  2001/10/21 19:04:55  hristov
+Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
+
 Revision 1.13  2001/07/23 12:01:30  hristov
 Initialisation added
 
@@ -66,8 +69,7 @@ Splitted from AliTPCtracking
 
 //_____________________________________________________________________________
 AliTPCtracker::AliTPCtracker(const AliTPCParam *par, Int_t eventn): 
-fkNIS(par->GetNInnerSector()/2), 
-fkNOS(par->GetNOuterSector()/2)
+AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
 {
   //---------------------------------------------------------------------
   // The main TPC tracker constructor
@@ -105,7 +107,7 @@ AliTPCtracker::~AliTPCtracker() {
   //------------------------------------------------------------------
   delete[] fInnerSec;
   delete[] fOuterSec;
-  delete fSeeds;
+  fSeeds->Delete(); delete fSeeds;
 }
 
 //_____________________________________________________________________________
@@ -481,7 +483,7 @@ void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
       for (Int_t js=0; js < nl+nm+nu; js++) {
        const AliTPCcluster *kcl;
         Double_t x2,   y2,   z2;
-        Double_t x3=0.,y3=0.;
+        Double_t x3=GetX(), y3=GetY(), z3=GetZ();
 
        if (js<nl) {
          const AliTPCRow& kr2=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
@@ -502,7 +504,7 @@ void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
             y2=xx2*sn+y2*cs;
          }
 
-        Double_t zz=z1 - z1/x1*(x1-x2); 
+        Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2); 
         if (TMath::Abs(zz-z2)>5.) continue;
 
         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
@@ -518,11 +520,11 @@ void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
        if (TMath::Abs(x[3]) > 1.2) continue;
        Double_t a=asin(x[2]);
        Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
-       if (TMath::Abs(zv)>10.) continue; 
+       if (TMath::Abs(zv-z3)>10.) continue; 
 
         Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
         Double_t sy2=kcl->GetSigmaY2(),     sz2=kcl->GetSigmaZ2();
-       Double_t sy3=100*0.025, sy=0.1, sz=0.1;
+       Double_t sy3=400*3./12., sy=0.1, sz=0.1;
 
        Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
        Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
@@ -618,7 +620,10 @@ Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
   }
 
   out->cd();
-  TTree tracktree("TPCf","Tree with TPC tracks");
+
+  char   tname[100];
+  sprintf(tname,"TreeT_TPC_%d",fEventN);
+  TTree tracktree(tname,"Tree with TPC tracks");
   AliTPCtrack *iotrack=0;
   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
 
@@ -645,7 +650,7 @@ Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
     }
     delete fSeeds->RemoveAt(i);
   }  
-  UnloadOuterSectors();
+  //UnloadOuterSectors();
 
   //tracking in inner sectors
   LoadInnerSectors();
@@ -666,6 +671,7 @@ Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
        if (FollowProlongation(t)) {
           if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
              t.CookdEdx();
+             CookLabel(pt,0.1); //For comparison only
              iotrack=pt;
              tracktree.Fill();
              UseClusters(&t,nc);
@@ -676,17 +682,9 @@ Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
     delete fSeeds->RemoveAt(i); 
   }  
   UnloadInnerSectors();
+  UnloadOuterSectors();
 
-  char   tname[100];
-  if (fEventN==-1) {
-    sprintf(tname,"TreeT_TPC");
-  }
-  else {
-    sprintf(tname,"TreeT_TPC_%d",fEventN);
-  }
-
-
-  tracktree.Write(tname);
+  tracktree.Write();
 
   cerr<<"Number of found tracks : "<<found<<endl;
 
@@ -983,6 +981,23 @@ void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
   for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
   dedx /= (nu-nl+1);
   SetdEdx(dedx);
+
+  //Very rough PID
+  Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
+
+  if (p<0.6) {
+    if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
+    if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
+    SetMass(0.93827); return;
+  }
+
+  if (p<1.2) {
+    if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
+    SetMass(0.93827); return;
+  }
+
+  SetMass(0.13957); return;
+
 }
 
 
index 04a3ff9..01a3edf 100644 (file)
@@ -135,9 +135,9 @@ private:
    Int_t fN;               //number of loaded sectors
    AliTPCSector *fSectors; //pointer to loaded sectors;
 
-  Int_t fEventN;                      //event number
+   Int_t fEventN;                      //event number
    AliTPCClustersArray fClustersArray; //array of TPC clusters
-   TObjArray *fSeeds;                  //array of track seeds 
+   TObjArray *fSeeds;                  //array of track seeds
 };
 
 #endif