From 4ab260d641203d4c5e08caa60f924b9650a2027f Mon Sep 17 00:00:00 2001
From: mariana <mariana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Tue, 19 Mar 2002 13:49:29 +0000
Subject: [PATCH] Last changes from I. Belikov - the tracking performs track
 refitting

---
 ITS/AliCascadeComparison.C |  10 +-
 ITS/AliCascadeVertexer.cxx |  10 ++
 ITS/AliITStrackerV2.cxx    | 187 ++++++++++++++++++++++++++++++-------
 ITS/AliITStrackerV2.h      |   1 +
 ITS/AliV0FindVertices.C    |   6 +-
 ITS/AliV0vertex.h          |  12 +++
 ITS/AliV0vertexer.cxx      |  20 +++-
 7 files changed, 202 insertions(+), 44 deletions(-)

diff --git a/ITS/AliCascadeComparison.C b/ITS/AliCascadeComparison.C
index 92758ff743f..bffb2938991 100644
--- a/ITS/AliCascadeComparison.C
+++ b/ITS/AliCascadeComparison.C
@@ -69,6 +69,7 @@ Int_t AliCascadeComparison(Int_t code=3312) {
        cTree->GetEvent(i);
        carray.AddLast(iovertex);
    }
+   delete cTree;
    cf->Close();
 
    /*** Check if the file with the "good" cascades exists ***/
@@ -129,7 +130,7 @@ Int_t AliCascadeComparison(Int_t code=3312) {
    hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
 
    Double_t mmin=cascadeMass-cascadeWindow, mmax=cascadeMass+cascadeWindow;
-   TH1F *cs =new TH1F("v0s","Cascade Effective Mass",40, mmin, mmax);
+   TH1F *cs =new TH1F("cs","Cascade Effective Mass",40, mmin, mmax);
    cs->SetXTitle("(GeV)"); cs->SetFillColor(6);
 
    Double_t pxg=0.,pyg=0.,ptg=0.;
@@ -264,7 +265,8 @@ Int_t AliCascadeComparison(Int_t code=3312) {
    c2->cd(2);
    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
    cs->SetXTitle("(GeV/c)"); 
-   cs->Fit("gaus","","",cascadeMass-cascadeWidth,cascadeMass+cascadeWidth);
+   //cs->Fit("gaus","","",cascadeMass-cascadeWidth,cascadeMass+cascadeWidth);
+   cs->Draw();
    Double_t max=cs->GetMaximum();
    TLine *line3 = 
       new TLine(cascadeMass-cascadeWidth,0.,cascadeMass-cascadeWidth,max);
@@ -357,7 +359,7 @@ Int_t good_cascades(GoodCascade *gc, Int_t max) {
       if (i==ngt) continue;
 
       /*** fiducial volume ***/
-      Double_t x=bp->Vx(), y=bp->Vy(), z=bp->Vz(), r2=x*x+y*y; //bachelor
+      Double_t x=bp->Vx(), y=bp->Vy(), r2=x*x+y*y; //bachelor
       if (r2<r2min) continue;
       if (r2>r2max) continue;
       TParticle *pp=gAlice->Particle(plab);
@@ -372,7 +374,7 @@ Int_t good_cascades(GoodCascade *gc, Int_t max) {
       gc[nc].code=code;
       gc[nc].plab=plab;   gc[nc].nlab=nlab; gc[nc].blab=blab;
       gc[nc].px=cp->Px(); gc[nc].py=cp->Py(); gc[nc].pz=cp->Pz();
-      gc[nc].x=x; gc[nc].y=y; gc[nc].z=z;
+      gc[nc].x=bp->Vx(); gc[nc].y=bp->Vy(); gc[nc].z=bp->Vz();
       nc++;
 
    }
diff --git a/ITS/AliCascadeVertexer.cxx b/ITS/AliCascadeVertexer.cxx
index c692598d831..980252316e3 100644
--- a/ITS/AliCascadeVertexer.cxx
+++ b/ITS/AliCascadeVertexer.cxx
@@ -80,6 +80,9 @@ AliCascadeVertexer::V0sTracks2CascadeVertices(const TFile *inp, TFile *out) {
        AliITStrackV2 *iotrack=new AliITStrackV2;
        branch->SetAddress(&iotrack);
        trkTree->GetEvent(i);
+
+       iotrack->PropagateTo(3.,0.0023,65.19); iotrack->PropagateTo(2.5,0.,0.);
+
        ntrack++; trks.AddLast(iotrack);
        
    }   
@@ -140,6 +143,13 @@ AliCascadeVertexer::V0sTracks2CascadeVertices(const TFile *inp, TFile *out) {
          if (r2 > fRmax*fRmax) continue;   // condition on fiducial zone
          if (r2 < fRmin*fRmin) continue;
 
+         {
+   //I.B.
+         Double_t x1,y1,z1; V0ver->GetXYZ(x1,y1,z1);
+         if (r2 > (x1*x1+y1*y1)) continue;
+         if (z*z > z1*z1) continue;
+         }
+
    // get cascade momentum
          
 	 Double_t px,py,pz; cascade.GetPxPyPz(px,py,pz);
diff --git a/ITS/AliITStrackerV2.cxx b/ITS/AliITStrackerV2.cxx
index 212af240653..591b95d2f1a 100644
--- a/ITS/AliITStrackerV2.cxx
+++ b/ITS/AliITStrackerV2.cxx
@@ -33,7 +33,7 @@
 //#define DEBUG
 
 #ifdef DEBUG
-Int_t LAB=5126;
+Int_t LAB=70201;
 #endif
 
 AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers
@@ -137,9 +137,9 @@ cout<<lay-1<<' '<<lad-1<<' '<<det-1<<' '<<c->GetY()<<' '<<c->GetZ()<<endl;
   fConstraint[0]=1; fConstraint[1]=0;
 }
 
-#ifdef DEBUG
+//#ifdef DEBUG
 static Int_t lbl;
-#endif
+//#endif
 
 Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) {
   //--------------------------------------------------------------------
@@ -178,7 +178,14 @@ Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) {
     nentr=(Int_t)tpcTree->GetEntries();
     for (Int_t i=0; i<nentr; i++) {
        tpcTree->GetEvent(i);
-       AliITStrackV2 *t=new AliITStrackV2(*itrack);
+       AliITStrackV2 *t=0;
+       try {
+           t=new AliITStrackV2(*itrack);
+       } catch (const Char_t *msg) {
+           cerr<<msg<<endl;
+           delete t;
+           continue;
+       }
        if (TMath::Abs(t->GetD())>4) continue;
 
        t->PropagateTo(80.,0.0053);
@@ -211,18 +218,14 @@ Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) {
        if (t==0) continue;           //this track has been already tracked
        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 
+lbl=tpcLabel;
 #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;
-       }
+       ResetTrackToFollow(*t);
        ResetBestTrack();
 
        for (FollowProlongation(); fI<kMaxLayer; fI++) {
@@ -233,14 +236,27 @@ cout<<tpcLabel<<" *****************\n";
 cout<<fBestTrack.GetNumberOfClusters()<<" number of clusters\n\n";
 #endif
 
-       if (fBestTrack.GetNumberOfClusters() >= kMaxLayer-kLayersToSkip) {
-          fBestTrack.SetLabel(tpcLabel);
-	  fBestTrack.CookdEdx();
-          CookLabel(&fBestTrack,0.); //For comparison only
-          itsTree.Fill();
-          UseClusters(&fBestTrack);
-          delete itsTracks.RemoveAt(i);
-       }
+      if (fBestTrack.GetNumberOfClusters() < kMaxLayer-kLayersToSkip) continue;
+
+        if (fConstraint[fPass]) {
+	   Int_t index[kMaxLayer];
+           Int_t k;
+           for (k=0; k<kMaxLayer; k++) index[k]=-1;
+           Int_t nc=fBestTrack.GetNumberOfClusters();
+           for (k=0; k<nc; k++) { 
+              Int_t idx=fBestTrack.GetClusterIndex(k),nl=(idx&0xf0000000)>>28;
+              index[nl]=idx; 
+           }
+           fBestTrack.~AliITStrackV2(); new(&fBestTrack) AliITStrackV2(*t);
+	   if (!RefitAt(3.7, &fBestTrack, index)) continue;
+	}
+
+        fBestTrack.SetLabel(tpcLabel);
+	fBestTrack.CookdEdx();
+        CookLabel(&fBestTrack,0.); //For comparison only
+        itsTree.Fill();
+        UseClusters(&fBestTrack);
+        delete itsTracks.RemoveAt(i);
 
      }
   }
@@ -534,10 +550,10 @@ Int_t AliITStrackerV2::TakeNextProlongation() {
   //  dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch 
   //--------------------------------------------------------------------
   AliITSlayer &layer=fLayers[fI];
-  AliITStrackV2 &t=fTracks[fI];
+  ResetTrackToFollow(fTracks[fI]);
 
-  Double_t dz=4*TMath::Sqrt(t.GetSigmaZ2() + kSigmaZ2[fI]);
-  Double_t dy=4*TMath::Sqrt(t.GetSigmaY2() + kSigmaY2[fI]);
+  Double_t dz=4*TMath::Sqrt(fTrackToFollow.GetSigmaZ2() + kSigmaZ2[fI]);
+  Double_t dy=4*TMath::Sqrt(fTrackToFollow.GetSigmaY2() + kSigmaY2[fI]);
 
   const AliITSclusterV2 *c=0; Int_t ci=-1;
   Double_t chi2=12345.;
@@ -545,15 +561,16 @@ Int_t AliITStrackerV2::TakeNextProlongation() {
     //if (c->GetLabel(0)!=TMath::Abs(lbl)) continue; 
     Int_t idet=c->GetDetectorIndex();
 
-    if (t.GetDetectorIndex()!=idet) {
+    if (fTrackToFollow.GetDetectorIndex()!=idet) {
        const AliITSdetector &det=layer.GetDetector(idet);
-       if (!t.Propagate(det.GetPhi(),det.GetR())) {
+       ResetTrackToFollow(fTracks[fI]);
+       if (!fTrackToFollow.Propagate(det.GetPhi(),det.GetR())) {
          //cerr<<"AliITStrackerV2::TakeNextProlongation: "
          //"propagation failed !\n";
          continue;
        }
-       t.SetDetectorIndex(idet);
-       if (TMath::Abs(t.GetZ()-GetZ()) > layer.GetR()+dz) continue;
+       fTrackToFollow.SetDetectorIndex(idet);
+       if (TMath::Abs(fTrackToFollow.GetZ()-GetZ())>layer.GetR()+dz) continue;
 
 #ifdef DEBUG
 cout<<fI<<" change detector !\n";
@@ -561,26 +578,29 @@ cout<<fI<<" change detector !\n";
 
     }
 
-    if (TMath::Abs(t.GetZ() - c->GetZ()) > dz) continue;
-    if (TMath::Abs(t.GetY() - c->GetY()) > dy) continue;
+    if (TMath::Abs(fTrackToFollow.GetZ() - c->GetZ()) > dz) continue;
+    if (TMath::Abs(fTrackToFollow.GetY() - c->GetY()) > dy) continue;
 
-    chi2=t.GetPredictedChi2(c); if (chi2<kMaxChi2) break;
+    chi2=fTrackToFollow.GetPredictedChi2(c); if (chi2<kMaxChi2) break;
   }
 
 #ifdef DEBUG
-cout<<fI<<" chi2="<<chi2<<' '<<t.GetY()<<' '<<t.GetZ()<<' '<<dy<<' '<<dz<<endl;
+cout<<fI<<" chi2="<<chi2<<' '
+    <<fTrackToFollow.GetY()<<' '<<fTrackToFollow.GetZ()<<' '
+    <<dy<<' '<<dz<<endl;
 #endif
 
   if (chi2>=kMaxChi2) return 0;
   if (!c) return 0;
 
-  ResetTrackToFollow(t);
-
   if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) {
      //cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
      return 0;
   }
 
+  if (fTrackToFollow.GetNumberOfClusters()>1)
+  if (TMath::Abs(fTrackToFollow.GetD())>4) return 0;
+
   fTrackToFollow.
     SetSampledEdx(c->GetQ(),fTrackToFollow.GetNumberOfClusters()-1); //b.b.
 
@@ -588,9 +608,10 @@ cout<<fI<<" chi2="<<chi2<<' '<<t.GetY()<<' '<<t.GetZ()<<' '<<dy<<' '<<dz<<endl;
    Double_t d=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ());
    fTrackToFollow.CorrectForMaterial(d);
   }
+
   if (fConstraint[fPass]) {
-     Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
-     fTrackToFollow.Improve(d,GetY(),GetZ());
+    Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
+    fTrackToFollow.Improve(d,GetY(),GetZ());
   }
 
 #ifdef DEBUG
@@ -831,5 +852,105 @@ Int_t AliITStrackerV2::AliITSlayer::InRoad() const {
   return ncl;
 }
 
+Bool_t AliITStrackerV2::RefitAt(Double_t x, AliITStrackV2 *t, Int_t *index) {
+  //--------------------------------------------------------------------
+  // This function refits a track at a given position
+  //--------------------------------------------------------------------
+  Int_t from, to, step;
+  if (x > t->GetX()) {
+      from=0; to=kMaxLayer;
+      step=+1;
+  } else {
+      from=kMaxLayer-1; to=-1;
+      step=-1;
+  }
+
+  for (Int_t i=from; i != to; i += step) {
+     AliITSlayer &layer=fLayers[i];
+     Double_t r=layer.GetR();
+ 
+     {
+     Double_t hI=i-0.5*step; 
+     if (hI==1.5 || hI==3.5) {             
+        Double_t rs=0.5*(fLayers[i-step].GetR() + r);
+        Double_t ds=0.011; if (hI==3.5) ds=0.0053;
+        if (!t->PropagateTo(rs,ds)) {
+          return kFALSE;
+        }
+     }
+     }
+
+     Double_t x,y,z;
+     if (!t->GetGlobalXYZat(r,x,y,z)) { 
+       return kFALSE;
+     }
+     Double_t phi=TMath::ATan2(y,x);
+     Int_t idet=layer.FindDetectorIndex(phi,z);
+     if (idet<0) { 
+       return kFALSE;
+     }
+     const AliITSdetector &det=layer.GetDetector(idet);
+     phi=det.GetPhi();
+     if (!t->Propagate(phi,det.GetR())) {
+       return kFALSE;
+     }
+     t->SetDetectorIndex(idet);
+
+     const AliITSclusterV2 *cl=0;
+     Double_t maxchi2=kMaxChi2;
+
+     Int_t idx=index[i];
+     if (idx>0) {
+        const AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(idx); 
+        if (idet != c->GetDetectorIndex()) {
+           idet=c->GetDetectorIndex();
+           const AliITSdetector &det=layer.GetDetector(idet);
+           if (!t->Propagate(det.GetPhi(),det.GetR())) {
+             return kFALSE;
+           }
+           t->SetDetectorIndex(idet);
+        }
+        Double_t chi2=t->GetPredictedChi2(c);
+        if (chi2<maxchi2) { cl=c; maxchi2=chi2; }
+        else return kFALSE;
+     }
+
+     /*
+     if (cl==0)
+     if (t->GetNumberOfClusters()>2) {
+        Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]);
+        Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
+        Double_t zmin=t->GetZ() - dz;
+        Double_t zmax=t->GetZ() + dz;
+        Double_t ymin=t->GetY() + phi*r - dy;
+        Double_t ymax=t->GetY() + phi*r + dy;
+        layer.SelectClusters(zmin,zmax,ymin,ymax);
+
+        const AliITSclusterV2 *c=0; Int_t ci=-1;
+        while ((c=layer.GetNextCluster(ci))!=0) {
+           if (idet != c->GetDetectorIndex()) continue;
+           Double_t chi2=t->GetPredictedChi2(c);
+           if (chi2<maxchi2) { cl=c; maxchi2=chi2; idx=ci; }
+        }
+     }
+     */
+
+     if (cl) {
+       if (!t->Update(cl,maxchi2,idx)) {
+          return kFALSE;
+       }
+     }
+
+     {
+     Double_t d=layer.GetThickness(t->GetY(),t->GetZ());
+     t->CorrectForMaterial(-step*d);
+     }
+
+  }
+
+  if (!t->PropagateTo(x,0.,0.)) return kFALSE;
+  return kTRUE;
+}
+
 
 
diff --git a/ITS/AliITStrackerV2.h b/ITS/AliITStrackerV2.h
index b6ff0392ae1..c0bca0a196f 100644
--- a/ITS/AliITStrackerV2.h
+++ b/ITS/AliITStrackerV2.h
@@ -27,6 +27,7 @@ public:
   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);
+  Bool_t RefitAt(Double_t x, AliITStrackV2 *t, Int_t *index);
 
   class AliITSdetector {
   public:
diff --git a/ITS/AliV0FindVertices.C b/ITS/AliV0FindVertices.C
index 6d1cf5c9f1e..b2e1cb57b16 100644
--- a/ITS/AliV0FindVertices.C
+++ b/ITS/AliV0FindVertices.C
@@ -15,10 +15,10 @@ Int_t AliV0FindVertices() {
    if (!in->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 2;}
 
    Double_t cuts[]={33,  // max. allowed chi2
-                    0.02,// min. allowed negative daughter's impact parameter 
-                    0.02,// min. allowed positive daughter's impact parameter 
+                    0.12,// min. allowed negative daughter's impact parameter 
+                    0.06,// min. allowed positive daughter's impact parameter 
                     0.080,// max. allowed DCA between the daughter tracks
-                    0.9985,// max. allowed cosine of V0's pointing angle
+                    0.998,// max. allowed cosine of V0's pointing angle
                     0.9,  // min. radius of the fiducial volume
                     2.9   // max. radius of the fiducial volume
                    };
diff --git a/ITS/AliV0vertex.h b/ITS/AliV0vertex.h
index 3de575a1331..2b6ebcdeeef 100644
--- a/ITS/AliV0vertex.h
+++ b/ITS/AliV0vertex.h
@@ -26,6 +26,8 @@ public:
   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 GetNPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const;
+  void GetPPxPyPz(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;}
@@ -49,6 +51,16 @@ private:
   ClassDef(AliV0vertex,1)   // reconstructed V0 vertex
 };
 
+inline 
+void AliV0vertex::GetNPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
+px=fNmom[0]; py=fNmom[1]; pz=fNmom[2];
+}
+
+inline 
+void AliV0vertex::GetPPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
+px=fPmom[0]; py=fPmom[1]; pz=fPmom[2];
+}
+
 #endif
 
 
diff --git a/ITS/AliV0vertexer.cxx b/ITS/AliV0vertexer.cxx
index 435d25ed379..cc8df89c1eb 100644
--- a/ITS/AliV0vertexer.cxx
+++ b/ITS/AliV0vertexer.cxx
@@ -64,6 +64,9 @@ Int_t AliV0vertexer::Tracks2V0vertices(const TFile *inp, TFile *out) {
        AliITStrackV2 *iotrack=new AliITStrackV2;
        branch->SetAddress(&iotrack);
        trkTree->GetEvent(i);
+
+       iotrack->PropagateTo(3.,0.0023,65.19); iotrack->PropagateTo(2.5,0.,0.);
+
        if (iotrack->Get1Pt() > 0.) {nneg++; negtrks.AddLast(iotrack);}
        else {npos++; postrks.AddLast(iotrack);}
    }   
@@ -79,11 +82,18 @@ Int_t AliV0vertexer::Tracks2V0vertices(const TFile *inp, TFile *out) {
       if (i%10==0) cerr<<nneg-i<<'\r';
       AliITStrackV2 *ntrk=(AliITStrackV2 *)negtrks.UncheckedAt(i);
 
+      if (TMath::Abs(ntrk->GetD())<fDPmin) continue;
+      if (TMath::Abs(ntrk->GetD())>fRmax) continue;
+
       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;
+         if (TMath::Abs(ptrk->GetD())>fRmax) continue;
+
+         if (TMath::Abs(ntrk->GetD())<fDNmin)
+         if (TMath::Abs(ptrk->GetD())<fDNmin) continue;
+
 
          AliITStrackV2 nt(*ntrk), pt(*ptrk), *pnt=&nt, *ppt=&pt;
 
@@ -106,7 +116,9 @@ Int_t AliV0vertexer::Tracks2V0vertices(const TFile *inp, TFile *out) {
          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;
+
+         if (cost < (5*fCPAmax-0.9-TMath::Sqrt(r2)*(fCPAmax-1))/4.1) continue;
+         //if (cost < fCPAmax) continue;
 
          //vertex.ChangeMassHypothesis(); //default is Lambda0 
 
@@ -270,8 +282,8 @@ Double_t AliV0vertexer::PropagateToDCA(AliITStrackV2 *n, AliITStrackV2 *p) {
     return 1.e+33;
   }  
 
-  //return TMath::Sqrt(dm);
-  return TMath::Sqrt(dx*dx + dy*dy + dz*dz);
+  return TMath::Sqrt(dm*TMath::Sqrt(dy2*dz2));
+  //return TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 }
 
 
-- 
2.43.0