]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliVertexerTracks.cxx
Replaced AliInfo with AliDebug
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracks.cxx
index 03cc3edcb956c67bb7982068f02f73e699d37eda..c88f2144ba39372f7c7f699c5941086d23c9c95b 100644 (file)
@@ -120,11 +120,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(const AliVEvent *vEvent)
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(AliVEvent *vEvent)
 {
 //
 // Primary vertex for current ESD or AOD event
@@ -176,11 +176,13 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
     }
     if(skipThis) continue;
 
+    // kITSrefit
+    if(fMode==0 && fITSrefit && !(track->GetStatus()&AliESDtrack::kITSrefit)) continue;
+
     if(!inputAOD) {  // ESD
       AliESDtrack* esdt = (AliESDtrack*)track;
       if(esdt->GetNcls(fMode) < fMinClusters) continue;
       if(fMode==0) {        // ITS mode
-       if(fITSrefit && !(esdt->GetStatus()&AliESDtrack::kITSrefit)) continue;
        Double_t x,p[5],cov[15];
        esdt->GetExternalParameters(x,p);
        esdt->GetExternalCovariance(cov);
@@ -206,7 +208,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const 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;
@@ -214,6 +216,17 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
     olddir->cd();
   }
 
+  // set vertex ID for tracks used in the fit
+  // (only for ESD)
+  if(!inputAOD) {
+    Int_t nIndices = fCurrentVertex->GetNIndices();
+    UShort_t *indices = fCurrentVertex->GetIndices();
+    for(Int_t ind=0; ind<nIndices; ind++) {
+      AliESDtrack *esdt = (AliESDtrack*)vEvent->GetTrack(indices[ind]);
+      esdt->SetVertexID(-1);
+    }
+  }
+
   return fCurrentVertex;
 }
 //----------------------------------------------------------------------------
@@ -248,7 +261,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
@@ -288,7 +301,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));
@@ -329,7 +342,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
       indices[jj] = fIdSel[jj];
     fCurrentVertex->SetIndices(nIndices,indices);
   }
-  delete [] indices; indices=NULL;
+  if (indices) {delete indices; indices=NULL;}
   //
 
   // set vertex title
@@ -344,9 +357,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;
@@ -580,6 +593,7 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
     } else { // neutral tracks (from a V0)
       track = new AliNeutralTrackParam(*(AliNeutralTrackParam*)trkArrayOrig.At(i));
     }
+
     // tgl cut
     if(TMath::Abs(track->GetTgl())>fMaxTgl) {
       AliDebug(1,Form(" rejecting track with tgl = %f",track->GetTgl()));
@@ -688,7 +702,7 @@ Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
     //
     Double_t ca=TMath::Cos(alphan-track->GetAlpha()), 
       sa=TMath::Sin(alphan-track->GetAlpha());
-    Double_t sf=track->GetSnp(), cf=TMath::Sqrt(1.- sf*sf);
+    Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
     Double_t sinNew =  sf*ca - cf*sa;
     if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
     if(!track->Rotate(alphan)) return kFALSE;
@@ -701,7 +715,7 @@ Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
 AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
                                                        TObjArray *trkArray,
                                                        UShort_t *id,
-                                                       Float_t *diamondxy) 
+                                                       Float_t *diamondxy) const
 {
 //
 // Removes tracks in trksTree from fit of inVtx
@@ -758,7 +772,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
     sumWi -= wWi;
     sumWiri -= wWiri;
 
-    // track chi2
+    // track contribution to chi2
     TMatrixD deltar = rv; deltar -= ri;
     TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
     Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
@@ -795,7 +809,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
   cov[5] = vVnew(2,2);
   
   // store data in the vertex object
-  AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks);
+  AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks+1); // the +1 is for the constraint
   outVtx->SetTitle(inVtx->GetTitle());
   UShort_t *inindices = inVtx->GetIndices();
   Int_t nIndices = nUsedTrks;
@@ -811,7 +825,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
     }
   }
   outVtx->SetIndices(nIndices,outindices);
-  delete [] outindices;
+  if (outindices) delete outindices;
 
   /*
     printf("Vertex before removing tracks:");
@@ -825,6 +839,106 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
   return outVtx;
 }
 //---------------------------------------------------------------------------
+AliESDVertex* AliVertexerTracks::RemoveConstraintFromVertex(AliESDVertex *inVtx,
+                                                    Float_t *diamondxyz,
+                                                    Float_t *diamondcov) const
+{
+//
+// Removes diamond constraint from fit of inVtx
+//
+
+  if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
+    printf("ERROR: primary vertex has no constraint: cannot remove it\n");
+    return 0x0;
+  }
+  if(inVtx->GetNContributors()<3) {
+    printf("ERROR: primary vertex has less than 2 tracks: cannot remove contraint\n");
+    return 0x0;
+  }
+
+  // diamond constraint 
+  TMatrixD vVb(3,3);
+  vVb(0,0) = diamondcov[0];
+  vVb(0,1) = diamondcov[1];
+  vVb(0,2) = 0.;
+  vVb(1,0) = diamondcov[1];
+  vVb(1,1) = diamondcov[2];
+  vVb(1,2) = 0.;
+  vVb(2,0) = 0.;
+  vVb(2,1) = 0.;
+  vVb(2,2) = diamondcov[5];
+  TMatrixD vVbInv(TMatrixD::kInverted,vVb);
+  TMatrixD rb(3,1);
+  rb(0,0) = diamondxyz[0];
+  rb(1,0) = diamondxyz[1];
+  rb(2,0) = diamondxyz[2];
+  TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
+
+  // input vertex
+  TMatrixD rv(3,1);
+  rv(0,0) = inVtx->GetXv();
+  rv(1,0) = inVtx->GetYv();
+  rv(2,0) = inVtx->GetZv();
+  TMatrixD vV(3,3);
+  Double_t cov[6];
+  inVtx->GetCovMatrix(cov);
+  vV(0,0) = cov[0];
+  vV(0,1) = cov[1]; vV(1,0) = cov[1];
+  vV(1,1) = cov[2];
+  vV(0,2) = cov[3]; vV(2,0) = cov[3];
+  vV(1,2) = cov[4]; vV(2,1) = cov[4]; 
+  vV(2,2) = cov[5];
+  TMatrixD vVInv(TMatrixD::kInverted,vV);
+  TMatrixD vVInvrv(vVInv,TMatrixD::kMult,rv);
+
+
+  TMatrixD sumWi = vVInv - vVbInv;
+
+
+  TMatrixD sumWiri = vVInvrv - vVbInvrb;
+
+  TMatrixD rvnew(3,1);
+  TMatrixD vVnew(3,3);
+
+  // new inverted of weights matrix
+  TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
+  vVnew = invsumWi;
+  // new position of primary vertex
+  rvnew.Mult(vVnew,sumWiri);
+
+  Double_t position[3];
+  position[0] = rvnew(0,0);
+  position[1] = rvnew(1,0);
+  position[2] = rvnew(2,0);
+  cov[0] = vVnew(0,0);
+  cov[1] = vVnew(0,1);
+  cov[2] = vVnew(1,1);
+  cov[3] = vVnew(0,2);
+  cov[4] = vVnew(1,2);
+  cov[5] = vVnew(2,2);
+
+
+  Double_t chi2 = inVtx->GetChi2();
+
+  // diamond constribution to chi2
+  TMatrixD deltar = rv; deltar -= rb;
+  TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
+  Double_t chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
+                   deltar(1,0)*vVbInvdeltar(1,0)+
+                   deltar(2,0)*vVbInvdeltar(2,0);
+  // remove from total chi2
+  chi2 -= chi2b;
+
+  // store data in the vertex object
+  AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,inVtx->GetNContributors()-1);
+  outVtx->SetTitle("VertexerTracksNoConstraint");
+  UShort_t *inindices = inVtx->GetIndices();
+  Int_t nIndices = inVtx->GetNIndices();
+  outVtx->SetIndices(nIndices,inindices);
+
+  return outVtx;
+}
+//---------------------------------------------------------------------------
 void AliVertexerTracks::SetCuts(Double_t *cuts) 
 {
 //
@@ -942,7 +1056,7 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
 {
   AliExternalTrackParam *track1;
   const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
-  static TClonesArray linarray("AliStrLine",knacc);
+  AliStrLine **linarray = new AliStrLine* [knacc];
   for(Int_t i=0; i<knacc; i++){
     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
     Double_t alpha=track1->GetAlpha();
@@ -955,7 +1069,6 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
     sigmasq[2]=track1->GetSigmaZ2();
     TMatrixD ri(3,1);
     TMatrixD wWi(3,3);
-    //if(!TrackToPoint(track1,ri,wWi)) {printf("WARNING\n");continue;}
     if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
     Double_t wmat[9];
     Int_t iel=0;
@@ -965,16 +1078,31 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
        iel++;
       }    
     }
-    new(linarray[i]) AliStrLine(pos,sigmasq,wmat,dir);
+    linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
   }
-  fVert=TrackletVertexFinder(&linarray,optUseWeights);
-  linarray.Clear("C");
+  fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
+  for(Int_t i=0; i<knacc; i++) delete linarray[i];
+  delete [] linarray;
 }
 //---------------------------------------------------------------------------
 AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
 {
-  // Calculate the point at minimum distance to prepared tracks 
+  // Calculate the point at minimum distance to prepared tracks (TClonesArray)
   const Int_t knacc = (Int_t)lines->GetEntriesFast();
+  AliStrLine** lines2 = new AliStrLine* [knacc];
+  for(Int_t i=0; i<knacc; i++){
+    lines2[i]= (AliStrLine*)lines->At(i);
+  }
+  AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights); 
+  delete [] lines2;
+  return vert;
+}
+
+//---------------------------------------------------------------------------
+AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
+{
+  // Calculate the point at minimum distance to prepared tracks (array of AliStrLine) 
+
   Double_t initPos[3]={0.,0.,0.};
 
   Double_t (*vectP0)[3]=new Double_t [knacc][3];
@@ -991,7 +1119,7 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t
   }
 
   for(Int_t i=0; i<knacc; i++){
-    AliStrLine* line1 = (AliStrLine*)lines->At(i)
+    AliStrLine *line1 = lines[i]
     Double_t p0[3],cd[3],sigmasq[3];
     Double_t wmat[9];
     if(!line1) printf("ERROR %d %d\n",i,knacc);
@@ -1072,6 +1200,7 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t
   delete vectP1;
   return theVert;
 }
+
 //---------------------------------------------------------------------------
 Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
                                       TMatrixD &ri,TMatrixD &wWi,
@@ -1209,8 +1338,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;
 }
@@ -1494,30 +1623,43 @@ void AliVertexerTracks::VertexFitter()
 AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(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");
 
@@ -1540,7 +1682,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
   Double_t d0z0[2],covd0z0[3];
   AliExternalTrackParam *t = 0;
   if(fCurrentVertex->GetNContributors()>0) {
-    indices = new UShort_t[fCurrentVertex->GetNContributors()];
+    indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
     for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
       indices[jj] = fIdSel[jj];
       t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
@@ -1557,16 +1699,18 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
   }
 
   // clean up
-  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
@@ -1581,9 +1725,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;
 }