]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliVertexerTracks.cxx
Bug fix
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracks.cxx
index 214dc4f2cb268d59f80983dd56204c0582f45742..1d1dae75e9fe5a759f1d6c5a8d03d4b46ee3c8ff 100644 (file)
@@ -42,6 +42,7 @@ AliVertexerTracks::AliVertexerTracks():
 TObject(),
 fVert(),
 fCurrentVertex(0),
+fFieldkG(-999.),
 fConstraint(kFALSE),
 fOnlyFitter(kFALSE),
 fMinTracks(1),
@@ -67,11 +68,12 @@ fDebug(0)
   SetNSigmad0(); 
   SetMaxd0z0(); 
 }
-//-----------------------------------------------------------------------------
-AliVertexerTracks::AliVertexerTracks(Double_t xStart, Double_t yStart):
+//----------------------------------------------------------------------------
+AliVertexerTracks::AliVertexerTracks(Double_t fieldkG):
 TObject(),
 fVert(),
 fCurrentVertex(0),
+fFieldkG(-999.),
 fConstraint(kFALSE),
 fOnlyFitter(kFALSE),
 fMinTracks(1),
@@ -90,12 +92,13 @@ fDebug(0)
 //
 // Standard constructor
 //
-  SetVtxStart(xStart,yStart);
+  SetVtxStart();
   SetVtxStartSigma();
   SetMinTracks();
   SetMinITSClusters();
   SetNSigmad0(); 
-  SetMaxd0z0(); 
+  SetMaxd0z0();
+  SetFieldkG(fieldkG);
 }
 //-----------------------------------------------------------------------------
 AliVertexerTracks::~AliVertexerTracks() 
@@ -119,9 +122,19 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
 //
   fCurrentVertex = 0;
 
+  // accept 1-track case only if constraint is available
+  if(!fConstraint && fMinTracks==1) fMinTracks=2;
+
+  // get Bz from ESD
+  //  SetFieldkG(esdEvent->GetMagneticField());
+
   // read tracks from ESD
   Int_t nTrksTot = (Int_t)esdEvent->GetNumberOfTracks();
-  if (nTrksTot<=0) return fCurrentVertex;
+  if (nTrksTot<=0) {
+    if(fDebug) printf("TooFewTracks\n");
+    TooFewTracks(esdEvent);
+    return fCurrentVertex;
+  } 
 
   TTree *trkTree = new TTree("TreeT","tracks");
   AliESDtrack *esdTrack = 0;
@@ -356,7 +369,7 @@ void AliVertexerTracks::OneTrackVertFinder()
   if(fDebug) printf("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArray.GetEntries());
   AliESDtrack *track1;
   track1 = (AliESDtrack*)fTrkArray.At(0);
-  Double_t field=GetField();
+  Double_t field=GetFieldkG();
   Double_t alpha=track1->GetAlpha();
   Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
   Double_t pos[3],dir[3]; 
@@ -389,7 +402,7 @@ void AliVertexerTracks::HelixVertexFinder()
   Double_t initPos[3];
   initPos[2] = 0.;
   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
-  Double_t field=GetField();
+  Double_t field=GetFieldkG();
 
   Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
 
@@ -472,7 +485,7 @@ Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree,Int_t optImpParCut)
   Double_t normdistx,normdisty;
   Float_t  d0z0[2],covd0z0[3]; 
   Double_t sigma;
-  Double_t field=GetField();
+  Double_t field=GetFieldkG();
 
   AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
 
@@ -494,9 +507,9 @@ Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree,Int_t optImpParCut)
     } else {              // optImpParCut==2
       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[1]);
+      normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
       if(normdistx < 3. && normdisty < 3. &&
-        (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[1]))) {
+        (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
        track->RelateToVertex(fCurrentVertex,field,100.);
       } else {
        track->RelateToVertex(initVertex,field,100.);
@@ -512,9 +525,9 @@ Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree,Int_t optImpParCut)
 
     if(fDebug) printf("trk %d; lab %d; |d0| = %f;  d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f\n",i,track->GetLabel(),TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),maxd0z0);
 
-    // during iteration 1, if fConstraint=kFALSE,
+    // during iterations 1 and 2, if fConstraint=kFALSE,
     // select tracks with d0oplusz0 < maxd0z0
-    if(optImpParCut==1 && !fConstraint && nEntries>=3 && 
+    if(optImpParCut>=1 && !fConstraint && nEntries>=3 && 
        fVert.GetNContributors()>0) {
       if(TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]) > maxd0z0) { 
        if(fDebug) printf("     rejected\n");
@@ -579,7 +592,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
     }
     Double_t alpha = track->GetAlpha();
     Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
-    track->AliExternalTrackParam::PropagateTo(xl,AliTracker::GetBz()); 
+    track->AliExternalTrackParam::PropagateTo(xl,GetFieldkG()); 
     // vector of track global coordinates
     TMatrixD ri(3,1);
     // covariance matrix of ri
@@ -687,17 +700,36 @@ void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
 }
 //---------------------------------------------------------------------------
 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
+{
+  AliESDtrack *track1;
+  Double_t field=GetFieldkG();
+  const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
+  TClonesArray *linarray = new TClonesArray("AliStrLine",1000);
+  TClonesArray &lines = *linarray;
+  for(Int_t i=0; i<knacc; i++){
+    track1 = (AliESDtrack*)fTrkArray.At(i);
+    Double_t alpha=track1->GetAlpha();
+    Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
+    Double_t pos[3],dir[3],sigmasq[3]; 
+    track1->GetXYZAt(mindist,field,pos);
+    track1->GetPxPyPzAt(mindist,field,dir);
+    sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
+    sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
+    sigmasq[2]=track1->GetSigmaZ2();
+    new(lines[i]) AliStrLine(pos,sigmasq,dir); 
+  }
+  fVert=TrackletVertexFinder(linarray,optUseWeights);
+  linarray->Delete();
+  delete linarray;
+}
+//---------------------------------------------------------------------------
+AliVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
 {
   // Calculate the point at minimum distance to prepared tracks 
   
-  Double_t initPos[3];
-  initPos[2] = 0.;
-  Double_t sigma=0;
-  for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
-  const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
-  Double_t field=GetField();
+  const Int_t knacc = (Int_t)lines->GetEntriesFast();
+  Double_t initPos[3]={0.,0.,0.};
 
-  AliESDtrack *track1;
   Double_t (*vectP0)[3]=new Double_t [knacc][3];
   Double_t (*vectP1)[3]=new Double_t [knacc][3];
   
@@ -706,19 +738,11 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
   for(Int_t i=0;i<3;i++)
     for(Int_t j=0;j<3;j++)sum[i][j]=0;
   for(Int_t i=0; i<knacc; i++){
-    track1 = (AliESDtrack*)fTrkArray.At(i);
-    Double_t alpha=track1->GetAlpha();
-    Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
-    Double_t pos[3],dir[3]; 
-    track1->GetXYZAt(mindist,field,pos);
-    track1->GetPxPyPzAt(mindist,field,dir);
-    AliStrLine *line1 = new AliStrLine(pos,dir); 
-    //    AliStrLine *line1 = new AliStrLine();
-    //    track1->ApproximateHelixWithLine(mindist,field,line1);
-
-    Double_t p0[3],cd[3];
+    AliStrLine* line1 = (AliStrLine*)lines->At(i); 
+    Double_t p0[3],cd[3],sigmasq[3];
     line1->GetP0(p0);
     line1->GetCd(cd);
+    line1->GetSigma2P0(sigmasq);
     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
     vectP0[i][0]=p0[0];
     vectP0[i][1]=p0[1];
@@ -730,23 +754,18 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
     Double_t matr[3][3];
     Double_t dknow[3];
     if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
-    if(optUseWeights==1){
-      Double_t sigmasq[3];
-      sigmasq[0]=track1->GetSigmaY2();
-      sigmasq[1]=track1->GetSigmaY2();
-      sigmasq[2]=track1->GetSigmaZ2();
-      GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
-    }
+    else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
+
 
     for(Int_t iii=0;iii<3;iii++){
       dsum[iii]+=dknow[iii]; 
       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
     }
-    delete line1;
   }
   
   Double_t vett[3][3];
   Double_t det=GetDeterminant3X3(sum);
+  Double_t sigma=0;
   
   if(det!=0){
     for(Int_t zz=0;zz<3;zz++){
@@ -769,14 +788,12 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
 
     sigma=TMath::Sqrt(sigma);
   }else{
-    Warning("StrLinVertexFinderMinDist","Finder did not succed");
     sigma=999;
   }
+  AliVertex theVert(initPos,sigma,knacc);
   delete vectP0;
   delete vectP1;
-  fVert.SetXYZ(initPos);
-  fVert.SetDispersion(sigma);
-  fVert.SetNContributors(knacc);
+  return theVert;
 }
 //---------------------------------------------------------------------------
 Bool_t AliVertexerTracks::TrackToPoint(AliESDtrack *t,
@@ -882,7 +899,7 @@ void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
   AliESDtrack *track2;
   Double_t pos[3],dir[3]; 
   Double_t alpha,mindist;
-  Double_t field=GetField();
+  Double_t field=GetFieldkG();
 
   for(Int_t i=0; i<nacc; i++){
     track1 = (AliESDtrack*)fTrkArray.At(i);
@@ -1052,7 +1069,7 @@ void AliVertexerTracks::VertexFitter(Bool_t useConstraint)
       alpha = t->GetAlpha();
       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
       // to vtxSeed (from finder)
-      t->AliExternalTrackParam::PropagateTo(xlStart,AliTracker::GetBz());   
+      t->AliExternalTrackParam::PropagateTo(xlStart,GetFieldkG());   
  
 
       // vector of track global coordinates
@@ -1123,6 +1140,8 @@ void AliVertexerTracks::VertexFitter(Bool_t useConstraint)
   covmatrix[4] = vV(1,2);
   covmatrix[5] = vV(2,2);
   
+  // for correct chi2/ndf, count constraint as additional "track"
+  if(fConstraint) nUsedTrks++;
   // store data in the vertex object
   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
 
@@ -1186,7 +1205,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TTree *trkTree,Bool_t o
       indices[jj] = (UShort_t)eta->GetID();
       if(optPropagate&&optUseFitter){
        if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
-         eta->RelateToVertex(fCurrentVertex,GetField(),100.);
+         eta->RelateToVertex(fCurrentVertex,GetFieldkG(),100.);
          if(fDebug) printf("Track %d propagated to found vertex\n",jj);
        }else{
          AliWarning("Found vertex outside beam pipe!");