When calculating a*a-b*b the form (a-b)*(a+b) is usually more numerically stable.
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracks.cxx
index de569ba..97a7626 100644 (file)
@@ -35,6 +35,9 @@
 #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,7 +70,8 @@ fITSrefit(kTRUE),
 fFiducialR(3.),
 fFiducialZ(30.),
 fnSigmaForUi00(1.5),
-fAlgo(1)
+fAlgo(1),
+fAlgoIter0(4)
 {
 //
 // Default constructor
@@ -100,14 +104,14 @@ fITSrefit(kTRUE),
 fFiducialR(3.),
 fFiducialZ(30.),
 fnSigmaForUi00(1.5),
-fAlgo(1)
+fAlgo(1),
+fAlgoIter0(4)
 {
 //
 // Standard constructor
 //
   SetVtxStart();
   SetVtxStartSigma();
-  SetTPCMode();
 }
 //-----------------------------------------------------------------------------
 AliVertexerTracks::~AliVertexerTracks() 
@@ -120,10 +124,10 @@ AliVertexerTracks::~AliVertexerTracks()
   if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
 }
 //----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
 {
 //
-// Primary vertex for current ESD event
+// Primary vertex for current ESD or AOD event
 // (Two iterations: 
 //  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
 //      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE  
@@ -131,11 +135,20 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
 //
   fCurrentVertex = 0;
 
+  TString evtype = vEvent->IsA()->GetName();
+  Bool_t inputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
+
+  if(inputAOD && fMode==1) {
+    printf("Error : AliVertexerTracks: no TPC-only vertex from AOD\n"); 
+    TooFewTracks(); 
+    return fCurrentVertex;
+  }
+
   // accept 1-track case only if constraint is available
   if(!fConstraint && fMinTracks==1) fMinTracks=2;
 
-  // read tracks from ESD
-  Int_t nTrks = (Int_t)esdEvent->GetNumberOfTracks();
+  // read tracks from AlivEvent
+  Int_t nTrks = (Int_t)vEvent->GetNumberOfTracks();
   if(nTrks<fMinTracks) {
     TooFewTracks();
     return fCurrentVertex;
@@ -150,38 +163,44 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
 
   Int_t nTrksOrig=0;
   AliExternalTrackParam *t=0;
-  // loop on ESD tracks
+  // loop on tracks
   for(Int_t i=0; i<nTrks; i++) {
-    AliESDtrack *esdt = esdEvent->GetTrack(i);
+    AliVTrack *track = (AliVTrack*)vEvent->GetTrack(i);
     // check tracks to skip
     Bool_t skipThis = kFALSE;
     for(Int_t j=0; j<fNTrksToSkip; j++) { 
-      if(esdt->GetID()==fTrksToSkip[j]) {
+      if(track->GetID()==fTrksToSkip[j]) {
        AliDebug(1,Form("skipping track: %d",i));
        skipThis = kTRUE;
       }
     }
     if(skipThis) continue;
 
-    // check number of clusters in ITS or TPC
-    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);
-      t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
-    } else if(fMode==1) { // TPC mode
-      t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
-      if(!t) continue;
-      Double_t radius = 2.8; //something less than the beam pipe radius
-      if(!PropagateTrackTo(t,radius)) 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);
+       t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
+      } else if(fMode==1) { // TPC mode
+       t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
+       if(!t) continue;
+       Double_t radius = 2.8; //something less than the beam pipe radius
+       if(!PropagateTrackTo(t,radius)) continue;
+      }
+    } else {          // AOD (only ITS mode)
+      Int_t ncls0=0;
+      for(Int_t l=0;l<6;l++) if(TESTBIT(track->GetITSClusterMap(),l)) ncls0++;
+      if(ncls0 < fMinClusters) continue;
+      t = new AliExternalTrackParam(track);
     }
     trkArrayOrig.AddLast(t);
-    idOrig[nTrksOrig]=(UShort_t)esdt->GetID();
+    idOrig[nTrksOrig]=(UShort_t)track->GetID();
     nTrksOrig++;
-  } // end loop on ESD tracks
+  } // end loop on tracks
   
   // call method that will reconstruct the vertex
   FindPrimaryVertex(&trkArrayOrig,idOrig);
@@ -232,7 +251,15 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
     if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
     Double_t cutsave = fDCAcut;  
     fDCAcut = fDCAcutIter0;
-    VertexFinder(1); // using weights, cutting dca < fDCAcutIter0
+    // vertex finder
+    switch (fAlgoIter0) {
+    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;  
+    }
     fDCAcut = cutsave;
     if(fVert.GetNContributors()>0) {
       fVert.GetXYZ(fNominalPos);
@@ -547,7 +574,12 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
 
   // loop on tracks
   for(Int_t i=0; i<nTrksOrig; i++) {
-    track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
+    AliExternalTrackParam *trackOrig=(AliExternalTrackParam*)trkArrayOrig.At(i);
+    if(trackOrig->Charge()!=0) { // normal tracks
+      track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
+    } 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()));
@@ -656,7 +688,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;
@@ -808,6 +840,8 @@ void AliVertexerTracks::SetCuts(Double_t *cuts)
   SetMinDetFitter(cuts[6]);
   SetMaxTgl(cuts[7]);
   SetFiducialRZ(cuts[8],cuts[9]);
+  fAlgo=(Int_t)(cuts[10]);
+  fAlgoIter0=(Int_t)(cuts[11]);
 
   return;
 }
@@ -821,7 +855,9 @@ void AliVertexerTracks::SetITSMode(Double_t dcacut,
                                   Double_t mindetfitter,
                                   Double_t maxtgl,
                                   Double_t fidR,
-                                  Double_t fidZ) 
+                                  Double_t fidZ,
+                                  Int_t finderAlgo,
+                                  Int_t finderAlgoIter0)
 {
 //
 //  Cut values for ITS mode
@@ -841,6 +877,8 @@ void AliVertexerTracks::SetITSMode(Double_t dcacut,
   SetMinDetFitter(mindetfitter);
   SetMaxTgl(maxtgl);
   SetFiducialRZ(fidR,fidZ);
+  fAlgo=finderAlgo;
+  fAlgoIter0=finderAlgoIter0;
 
   return; 
 }
@@ -854,7 +892,9 @@ void AliVertexerTracks::SetTPCMode(Double_t dcacut,
                                   Double_t mindetfitter,
                                   Double_t maxtgl,
                                   Double_t fidR,
-                                  Double_t fidZ) 
+                                  Double_t fidZ,
+                                  Int_t finderAlgo,
+                                  Int_t finderAlgoIter0) 
 {
 //
 //  Cut values for TPC mode
@@ -870,6 +910,8 @@ void AliVertexerTracks::SetTPCMode(Double_t dcacut,
   SetMinDetFitter(mindetfitter);
   SetMaxTgl(maxtgl);
   SetFiducialRZ(fidR,fidZ);
+  fAlgo=finderAlgo;
+  fAlgoIter0=finderAlgoIter0;
 
   return; 
 }
@@ -913,7 +955,8 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
     sigmasq[2]=track1->GetSigmaZ2();
     TMatrixD ri(3,1);
     TMatrixD wWi(3,3);
-    if(!TrackToPoint(track1,ri,wWi)) continue;
+    //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;
     for(Int_t ia=0;ia<3;ia++){
@@ -922,7 +965,7 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
        iel++;
       }    
     }
-    new(linarray[i]) AliStrLine(pos,sigmasq,wmat,dir);     
+    new(linarray[i]) AliStrLine(pos,sigmasq,wmat,dir);
   }
   fVert=TrackletVertexFinder(&linarray,optUseWeights);
   linarray.Clear("C");
@@ -931,7 +974,6 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
 AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
 {
   // Calculate the point at minimum distance to prepared tracks 
-  
   const Int_t knacc = (Int_t)lines->GetEntriesFast();
   Double_t initPos[3]={0.,0.,0.};
 
@@ -952,6 +994,7 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t
     AliStrLine* line1 = (AliStrLine*)lines->At(i); 
     Double_t p0[3],cd[3],sigmasq[3];
     Double_t wmat[9];
+    if(!line1) printf("ERROR %d %d\n",i,knacc);
     line1->GetP0(p0);
     line1->GetCd(cd);
     line1->GetSigma2P0(sigmasq);