]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliVertexerTracks.cxx
Init the TList with refence objects when using the copy ctor, added protection agains...
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracks.cxx
index fc87a24f85b514875bc11caa695ea625a8c45db1..5a68582d75cb0818ac5aab79f485ba3476b5e26c 100644 (file)
 
 //---- Root headers --------
 #include <TSystem.h>
+#include <TClonesArray.h>
 #include <TDirectory.h>
 #include <TFile.h>
-#include <TTree.h>
-#include <TMatrixD.h>
 //---- AliRoot headers -----
-#include "AliLog.h"
 #include "AliStrLine.h"
 #include "AliExternalTrackParam.h"
 #include "AliNeutralTrackParam.h"
 #include "AliVEvent.h"
 #include "AliVTrack.h"
-#include "AliESDEvent.h"
 #include "AliESDtrack.h"
 #include "AliVertexerTracks.h"
 
@@ -67,6 +64,7 @@ fMaxd0z0(0.5),
 fMinDetFitter(100.),
 fMaxTgl(1000.),
 fITSrefit(kTRUE),
+fITSpureSA(kFALSE),
 fFiducialR(3.),
 fFiducialZ(30.),
 fnSigmaForUi00(1.5),
@@ -101,6 +99,7 @@ fMaxd0z0(0.5),
 fMinDetFitter(100.),
 fMaxTgl(1000.),
 fITSrefit(kTRUE),
+fITSpureSA(kFALSE),
 fFiducialR(3.),
 fFiducialZ(30.),
 fnSigmaForUi00(1.5),
@@ -120,11 +119,11 @@ AliVertexerTracks::~AliVertexerTracks()
   // The objects pointed by the following pointer are not owned
   // by this class and are not deleted
   fCurrentVertex = 0;
-  if(fTrksToSkip) { delete fTrksToSkip; fTrksToSkip=NULL; }
-  if(fIdSel) { delete fIdSel; fIdSel=NULL; }
+  if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
+  if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
 }
 //----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::FindPrimaryVertex(AliVEvent *vEvent)
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
 {
 //
 // Primary vertex for current ESD or AOD event
@@ -176,6 +175,11 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(AliVEvent *vEvent)
     }
     if(skipThis) continue;
 
+    // skip pure ITS SA tracks (default)
+    if(!fITSpureSA && (track->GetStatus()&AliESDtrack::kITSpureSA)) continue;
+    // or use only pure ITS SA tracks
+    if(fITSpureSA && !(track->GetStatus()&AliESDtrack::kITSpureSA)) continue;
+
     // kITSrefit
     if(fMode==0 && fITSrefit && !(track->GetStatus()&AliESDtrack::kITSrefit)) continue;
 
@@ -208,7 +212,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(AliVEvent *vEvent)
   FindPrimaryVertex(&trkArrayOrig,idOrig);
 
   if(fMode==0) trkArrayOrig.Delete();
-  delete idOrig; idOrig=NULL;
+  delete [] idOrig; idOrig=NULL;
 
   if(f) {
     f->Close(); delete f; f = NULL;
@@ -230,7 +234,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(AliVEvent *vEvent)
   return fCurrentVertex;
 }
 //----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const TObjArray *trkArrayOrig,
                                                   UShort_t *idOrig)
 {
 //
@@ -261,7 +265,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
     // fill fTrkArraySel, for VertexFinder()
     fIdSel = new UShort_t[nTrksOrig];
     PrepareTracks(*trkArrayOrig,idOrig,0);
-    if(fIdSel) { delete fIdSel; fIdSel=NULL; }
+    if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
     Double_t cutsave = fDCAcut;  
     fDCAcut = fDCAcutIter0;
     // vertex finder
@@ -301,7 +305,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
   //                   between initVertex and fCurrentVertex) 
   for(Int_t iter=1; iter<=2; iter++) {
     if(fOnlyFitter && iter==1) continue; 
-    if(fIdSel) { delete fIdSel; fIdSel=NULL; }
+    if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
     fIdSel = new UShort_t[nTrksOrig];
     Int_t nTrksSel = PrepareTracks(*trkArrayOrig,idOrig,iter);
     AliDebug(1,Form("N tracks selected in iteration %d: %d",iter,nTrksSel));
@@ -342,7 +346,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
       indices[jj] = fIdSel[jj];
     fCurrentVertex->SetIndices(nIndices,indices);
   }
-  if (indices) {delete indices; indices=NULL;}
+  if (indices) {delete [] indices; indices=NULL;}
   //
 
   // set vertex title
@@ -357,9 +361,9 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
   AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
 
   // clean up
-  delete fIdSel; fIdSel=NULL;
+  delete [] fIdSel; fIdSel=NULL;
   fTrkArraySel.Delete();
-  if(fTrksToSkip) { delete fTrksToSkip; fTrksToSkip=NULL; }
+  if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
   //
   
   return fCurrentVertex;
@@ -372,7 +376,7 @@ Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3])
  return det;
 }
 //-------------------------------------------------------------------------
-void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d)
+void AliVertexerTracks::GetStrLinDerivMatrix(const Double_t *p0,const Double_t *p1,Double_t (*m)[3],Double_t *d)
 {
   //
   Double_t x12=p0[0]-p1[0];
@@ -394,7 +398,7 @@ void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t
 
 }
 //--------------------------------------------------------------------------  
-void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
+void AliVertexerTracks::GetStrLinDerivMatrix(const Double_t *p0,const Double_t *p1,const Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
 {
   //
   Double_t x12=p1[0]-p0[0];
@@ -437,7 +441,7 @@ void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t
 
   }
 //--------------------------------------------------------------------------   
-Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0)
+Double_t AliVertexerTracks::GetStrLinMinDist(const Double_t *p0,const Double_t *p1,const Double_t *x0)
 {
   //
   Double_t x12=p0[0]-p1[0];
@@ -562,8 +566,8 @@ void AliVertexerTracks::HelixVertexFinder()
   fVert.SetNContributors(ncombi);
 }
 //----------------------------------------------------------------------------
-Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
-                                      UShort_t *idOrig,
+Int_t AliVertexerTracks::PrepareTracks(const TObjArray &trkArrayOrig,
+                                      const UShort_t *idOrig,
                                       Int_t optImpParCut) 
 {
 //
@@ -600,7 +604,7 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
       delete track; continue;
     }
 
-    Bool_t propagateOK = kFALSE;
+    Bool_t propagateOK = kFALSE, cutond0z0 = kTRUE;
     // propagate track to vertex
     if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
       propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
@@ -608,11 +612,15 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
       fCurrentVertex->GetSigmaXYZ(sigmaCurr);
       normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
       normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
+      AliDebug(1,Form("normdistx %f  %f    %f",fCurrentVertex->GetXv(),fNominalPos[0],TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0])));
+      AliDebug(1,Form("normdisty %f  %f    %f",fCurrentVertex->GetYv(),fNominalPos[1],TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2])));
+      AliDebug(1,Form("sigmaCurr %f %f    %f",sigmaCurr[0],sigmaCurr[1],TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2])));
       if(normdistx < 3. && normdisty < 3. &&
         (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
        propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
       } else {
        propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
+       if(fConstraint) cutond0z0=kFALSE;
       }
     }
 
@@ -625,8 +633,7 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
     maxd0rphi = fNSigma*sigmad0;
     if(optImpParCut==1) maxd0rphi *= 5.;
     maxd0rphi = TMath::Min(maxd0rphi,fFiducialR); 
-    //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]); // for future improvement
-    //maxd0z0 = 10.*fNSigma*sigmad0z0;
+    //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]);
 
     AliDebug(1,Form("trk %d; id %d; |d0| = %f;  d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f",i,(Int_t)idOrig[i],TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),fMaxd0z0));
 
@@ -649,7 +656,7 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
     // if fConstraint=kTRUE, during iteration 2,
     // select tracks with d0oplusz0 < fMaxd0z0
     if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
-       ( fConstraint && optImpParCut==2)) {
+       ( fConstraint && optImpParCut==2 && cutond0z0)) {
       if(nTrksOrig>=3 && 
         TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) { 
        AliDebug(1,"     rejected");
@@ -713,9 +720,9 @@ Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
 }
 //---------------------------------------------------------------------------
 AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
-                                                       TObjArray *trkArray,
+                                                       const TObjArray *trkArray,
                                                        UShort_t *id,
-                                                       Float_t *diamondxy) const
+                                                       const Float_t *diamondxy) const
 {
 //
 // Removes tracks in trksTree from fit of inVtx
@@ -825,7 +832,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
     }
   }
   outVtx->SetIndices(nIndices,outindices);
-  if (outindices) delete outindices;
+  if (outindices) delete [] outindices;
 
   /*
     printf("Vertex before removing tracks:");
@@ -1030,7 +1037,7 @@ void AliVertexerTracks::SetTPCMode(Double_t dcacut,
   return; 
 }
 //---------------------------------------------------------------------------
-void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) 
+void AliVertexerTracks::SetSkipTracks(Int_t n,const Int_t *skipped) 
 {
 //
 // Mark the tracks not to be used in the vertex reconstruction.
@@ -1085,7 +1092,7 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
   delete [] linarray;
 }
 //---------------------------------------------------------------------------
-AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
+AliESDVertex AliVertexerTracks::TrackletVertexFinder(const TClonesArray *lines, Int_t optUseWeights)
 {
   // Calculate the point at minimum distance to prepared tracks (TClonesArray)
   const Int_t knacc = (Int_t)lines->GetEntriesFast();
@@ -1171,7 +1178,7 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const I
   Double_t det=GetDeterminant3X3(sum);
   Double_t sigma=0;
   
-  if(det!=0){
+  if(TMath::Abs(det) > kAlmost0){
     for(Int_t zz=0;zz<3;zz++){
       for(Int_t ww=0;ww<3;ww++){
        for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
@@ -1196,8 +1203,8 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const I
   }
   AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
   theVert.SetDispersion(sigma);
-  delete vectP0;
-  delete vectP1;
+  delete [] vectP0;
+  delete [] vectP1;
   return theVert;
 }
 
@@ -1338,8 +1345,8 @@ void AliVertexerTracks::TooFewTracks()
   }
 
   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete(); 
-  if(fIdSel) {delete fIdSel; fIdSel=NULL;}
-  if(fTrksToSkip) {delete fTrksToSkip; fTrksToSkip=NULL;}
+  if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
+  if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
 
   return;
 }
@@ -1620,33 +1627,46 @@ void AliVertexerTracks::VertexFitter()
   return;
 }
 //----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
+AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArray,
                                                         UShort_t *id,
                                                         Bool_t optUseFitter,
-                                                        Bool_t optPropagate) 
+                                                        Bool_t optPropagate,
+                                                        Bool_t optUseDiamondConstraint) 
 {
 //
 // Return vertex from tracks (AliExternalTrackParam) in array
 //
   fCurrentVertex = 0;
-  SetConstraintOff();
+  // set optUseDiamondConstraint=TRUE only if you are reconstructing the 
+  // primary vertex!
+  if(optUseDiamondConstraint) {
+    SetConstraintOn();
+  } else {    
+    SetConstraintOff();
+  }
 
   // get tracks and propagate them to initial vertex position
   fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
   Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
-  if(nTrksSel <  TMath::Max(2,fMinTracks)) {
+  if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) ||
+     (optUseDiamondConstraint && nTrksSel<1)) {
     TooFewTracks();
     return fCurrentVertex;
   }
  
   // vertex finder
-  switch (fAlgo) {
-    case 1: StrLinVertexFinderMinDist(1); break;
-    case 2: StrLinVertexFinderMinDist(0); break;
-    case 3: HelixVertexFinder();          break;
-    case 4: VertexFinder(1);              break;
-    case 5: VertexFinder(0);              break;
-    default: printf("Wrong algorithm\n"); break;  
+  if(nTrksSel==1) {
+    AliDebug(1,"Just one track");
+    OneTrackVertFinder();
+  } else {
+    switch (fAlgo) {
+      case 1: StrLinVertexFinderMinDist(1); break;
+      case 2: StrLinVertexFinderMinDist(0); break;
+      case 3: HelixVertexFinder();          break;
+      case 4: VertexFinder(1);              break;
+      case 5: VertexFinder(0);              break;
+      default: printf("Wrong algorithm\n"); break;  
+    }
   }
   AliDebug(1," Vertex finding completed\n");
 
@@ -1686,16 +1706,18 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
   }
 
   // clean up
-  if (indices) {delete indices; indices=NULL;}
-  delete fIdSel; fIdSel=NULL;
+  if (indices) {delete [] indices; indices=NULL;}
+  delete [] fIdSel; fIdSel=NULL;
   fTrkArraySel.Delete();
   
   return fCurrentVertex;
 }
 //----------------------------------------------------------------------------
 AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
-                                                        Bool_t optUseFitter,
-                                                        Bool_t optPropagate) 
+                                                           Bool_t optUseFitter,
+                                                           Bool_t optPropagate,
+                                                           Bool_t optUseDiamondConstraint)
+
 {
 //
 // Return vertex from array of ESD tracks
@@ -1710,9 +1732,9 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
     id[i] = (UShort_t)esdt->GetID();
   }
     
-  VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate);
+  VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint);
 
-  delete id; id=NULL;
+  delete [] id; id=NULL;
 
   return fCurrentVertex;
 }