]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/ESD/AliVertexerTracks.cxx
In AddTaskPHOSPi0Flow.C set Cent. Bin past event buffers/lists,
[u/mrichter/AliRoot.git] / STEER / ESD / AliVertexerTracks.cxx
index 6e482c9f02948c83ee4a6341fee26692690bc844..b5619809a1ae89ff2fd7274aba975590b8d8cc24 100644 (file)
@@ -58,7 +58,7 @@ fNTrksToSkip(0),
 fConstraint(kFALSE),
 fOnlyFitter(kFALSE),
 fMinTracks(1),
-fMinClusters(5),
+fMinClusters(3),
 fDCAcut(0.1),
 fDCAcutIter0(0.1),
 fNSigma(3.),
@@ -86,7 +86,10 @@ fMVScanStep(3.),
 fMVMaxWghNtr(10.),
 fMVFinalWBinary(kTRUE),
 fBCSpacing(50),
-fMVVertices(0)
+fMVVertices(0),
+fClusterize(kFALSE),
+fDeltaZCutForCluster(0.1),
+fnSigmaZCutForCluster(999999.)
 {
 //
 // Default constructor
@@ -108,7 +111,7 @@ fNTrksToSkip(0),
 fConstraint(kFALSE),
 fOnlyFitter(kFALSE),
 fMinTracks(1),
-fMinClusters(5),
+fMinClusters(3),
 fDCAcut(0.1),
 fDCAcutIter0(0.1),
 fNSigma(3.),
@@ -136,7 +139,10 @@ fMVScanStep(3.),
 fMVMaxWghNtr(10.),
 fMVFinalWBinary(kTRUE),
 fBCSpacing(50),
-fMVVertices(0)
+fMVVertices(0),
+fClusterize(kFALSE),
+fDeltaZCutForCluster(0.1),
+fnSigmaZCutForCluster(999999.)
 {
 //
 // Standard constructor
@@ -192,6 +198,8 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
   if(nTrks>500) f = new TFile("VertexerTracks.root","recreate");
   TObjArray trkArrayOrig(nTrks);
   UShort_t *idOrig = new UShort_t[nTrks];
+  Double_t *zTr = new Double_t[nTrks];
+  Double_t *err2zTr = new Double_t[nTrks];
 
   Int_t nTrksOrig=0;
   AliExternalTrackParam *t=0;
@@ -234,7 +242,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
       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);
+      t = new AliExternalTrackParam(); t->CopyFromVTrack(track);
     }
 
     // use TOF info about bunch crossing
@@ -246,16 +254,21 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
     //
     trkArrayOrig.AddLast(t);
     idOrig[nTrksOrig]=(UShort_t)track->GetID();
+    zTr[nTrksOrig]=t->GetZ();
+    err2zTr[nTrksOrig]=t->GetSigmaZ2();
+
     nTrksOrig++;
   } // end loop on tracks
   
   // call method that will reconstruct the vertex
-  FindPrimaryVertex(&trkArrayOrig,idOrig);
-
+  if(fClusterize) FindAllVertices(nTrksOrig,&trkArrayOrig,zTr,err2zTr,idOrig);
+  else FindPrimaryVertex(&trkArrayOrig,idOrig);
   if(!inputAOD) AnalyzePileUp((AliESDEvent*)vEvent);
 
   if(fMode==0) trkArrayOrig.Delete();
   delete [] idOrig; idOrig=NULL;
+  delete [] zTr; zTr=NULL;
+  delete [] err2zTr; err2zTr=NULL;
 
   if(f) {
     f->Close(); delete f; f = NULL;
@@ -273,7 +286,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
       esdt->SetVertexID(-1);
     }
   }
-
   return fCurrentVertex;
 }
 //----------------------------------------------------------------------------
@@ -317,7 +330,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const TObjArray *trkArrayOrig
     case 3: HelixVertexFinder();          break;
     case 4: VertexFinder(1);              break;
     case 5: VertexFinder(0);              break;
-    default: printf("Wrong algorithm\n"); break;  
+    default: {AliFatal(Form("Wrong seeder algorithm %d",fAlgoIter0));} break;  
     }
     fDCAcut = cutsave;
     if(fVert.GetNContributors()>0) {
@@ -371,7 +384,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const TObjArray *trkArrayOrig
         case kVertexFinder1             : VertexFinder(1);              break;
         case kVertexFinder0             : VertexFinder(0);              break;
        case kMultiVertexer             : FindVerticesMV(); multiMode = kTRUE; break;
-        default: printf("Wrong algorithm"); break;  
+        default: {AliFatal(Form("Wrong vertexer algorithm %d",fAlgo));} break;  
        }
       }
       AliDebug(1," Vertex finding completed");
@@ -433,7 +446,7 @@ void AliVertexerTracks::GetStrLinDerivMatrix(const Double_t *p0,const Double_t *
   m[1][1]=2-2/kk*y12*y12;
   m[1][2]=-2/kk*y12*z12;
   m[2][0]=-2/kk*x12*z12;
-  m[2][1]=-2*y12*z12;
+  m[2][1]=-2/kk*y12*z12;
   m[2][2]=2-2/kk*z12*z12;
   d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
   d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
@@ -868,7 +881,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
   for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
     Bool_t copyindex=kTRUE;
     for(Int_t l=0; l<ntrks; l++) {
-      if(inindices[k]==id[l]) copyindex=kFALSE;
+      if(inindices[k]==id[l]) {copyindex=kFALSE; break;}
     }
     if(copyindex) {
       outindices[j] = inindices[k]; j++;
@@ -998,11 +1011,11 @@ void AliVertexerTracks::SetCuts(Double_t *cuts, Int_t ncuts)
   if (ncuts>1) SetDCAcutIter0(cuts[1]);
   if (ncuts>2) SetMaxd0z0(cuts[2]);
   if (ncuts>3) if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
-  if (ncuts>4) SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
-  if (ncuts>5) SetMinTracks((Int_t)(cuts[4]));
-  if (ncuts>6) SetNSigmad0(cuts[5]);
-  if (ncuts>7) SetMinDetFitter(cuts[6]);
-  if (ncuts>8) SetMaxTgl(cuts[7]);
+  if (ncuts>3) SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
+  if (ncuts>4) SetMinTracks((Int_t)(cuts[4]));
+  if (ncuts>5) SetNSigmad0(cuts[5]);
+  if (ncuts>6) SetMinDetFitter(cuts[6]);
+  if (ncuts>7) SetMaxTgl(cuts[7]);
   if (ncuts>9) SetFiducialRZ(cuts[8],cuts[9]);
   if (ncuts>10) fAlgo=(Int_t)(cuts[10]);
   if (ncuts>11) fAlgoIter0=(Int_t)(cuts[11]);
@@ -1017,10 +1030,13 @@ void AliVertexerTracks::SetCuts(Double_t *cuts, Int_t ncuts)
   if (ncuts>19) SetMVFinalWBinary(cuts[19]>0);
   if (ncuts>20) if (cuts[20]>20.)  SetBCSpacing(int(cuts[20]));
   //
-  if (ncuts>21) {
-    if (fAlgo==kMultiVertexer) SetSelectOnTOFBunchCrossing(kTRUE,kTRUE);
-    else                       SetSelectOnTOFBunchCrossing(kFALSE,kTRUE);
-  }
+  if (ncuts>21) if (cuts[21]>0.5)  SetUseTrackClusterization(kTRUE);
+  if (ncuts>22) SetDeltaZCutForCluster(cuts[22]);
+  if (ncuts>23) SetnSigmaZCutForCluster(cuts[23]);  
+  //
+  if (fAlgo==kMultiVertexer || fClusterize) SetSelectOnTOFBunchCrossing(kTRUE,kTRUE);
+  else                       SetSelectOnTOFBunchCrossing(kFALSE,kTRUE);
+  //
   return;
 }
 //---------------------------------------------------------------------------
@@ -1133,7 +1149,10 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
     sigmasq[2]=track1->GetSigmaZ2();
     TMatrixD ri(3,1);
     TMatrixD wWi(3,3);
-    if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
+    if(!TrackToPoint(track1,ri,wWi)) {
+      optUseWeights=kFALSE;
+      AliDebug(1,"WARNING: TrackToPoint failed");
+    }
     Double_t wmat[9];
     Int_t iel=0;
     for(Int_t ia=0;ia<3;ia++){
@@ -1186,7 +1205,7 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const 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);
+    if(!line1) {printf("ERROR %d %d\n",i,knacc); continue;}
     line1->GetP0(p0);
     line1->GetCd(cd);
     line1->GetSigma2P0(sigmasq);
@@ -1778,6 +1797,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArr
     }
   }
   AliDebug(1," Vertex finding completed\n");
+  Double_t vdispersion=fVert.GetDispersion();
 
   // vertex fitter
   if(optUseFitter) {
@@ -1790,7 +1810,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArr
     Int_t    nTrksUsed=fVert.GetNContributors();
     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);    
   }
-  fCurrentVertex->SetDispersion(fVert.GetDispersion());
+  fCurrentVertex->SetDispersion(vdispersion);
 
 
   // set indices of used tracks and propagate track to found vertex
@@ -1938,7 +1958,7 @@ void AliVertexerTracks::FindVerticesMV()
   double step = fMVScanStep>1 ?  fMVScanStep : 1.;
   double zmx = 3*TMath::Sqrt(fNominalCov[5]);
   double zmn = -zmx;
-  int nz = (int)(zmx-zmn)/step; if (nz<1) nz=1;
+  int nz = TMath::Nint((zmx-zmn)/step); if (nz<1) nz=1;
   double dz = (zmx-zmn)/nz;
   int izStart=0;
   AliDebug(2,Form("%d seeds between %f and %f",nz,zmn+dz/2,zmx+dz/2));
@@ -2072,3 +2092,113 @@ void AliVertexerTracks::AnalyzePileUp(AliESDEvent* esdEv)
   //
 }
 
+//______________________________________________________
+void AliVertexerTracks::FindAllVertices(Int_t nTrksOrig, 
+                                       const TObjArray *trkArrayOrig,
+                                       Double_t* zTr, 
+                                       Double_t* err2zTr, 
+                                       UShort_t* idOrig){
+
+  // clusterize tracks using z coordinates of intersection with beam axis
+  // and compute all vertices 
+
+  UShort_t* posOrig=new UShort_t[nTrksOrig];
+  for(Int_t iTr=0; iTr<nTrksOrig; iTr++) posOrig[iTr]=iTr;
+
+  // sort points along Z
+  AliDebug(1,Form("Sort points along Z, used tracks %d",nTrksOrig));
+  for(Int_t iTr1=0; iTr1<nTrksOrig; iTr1++){
+    for(Int_t iTr2=iTr1+1; iTr2<nTrksOrig; iTr2++){
+      if(zTr[iTr1]>zTr[iTr2]){
+       Double_t tmpz=zTr[iTr2];
+       Double_t tmperr=err2zTr[iTr2];
+       UShort_t tmppos=posOrig[iTr2];
+       UShort_t tmpid=idOrig[iTr2];
+       zTr[iTr2]=zTr[iTr1];
+       err2zTr[iTr2]=err2zTr[iTr1];
+       posOrig[iTr2]=posOrig[iTr1];
+       idOrig[iTr2]=idOrig[iTr1];
+       zTr[iTr1]=tmpz;
+       err2zTr[iTr1]=tmperr;
+       idOrig[iTr1]=tmpid;
+       posOrig[iTr1]=tmppos;
+      }
+    }
+  }
+
+  // clusterize
+  Int_t nClusters=0;
+  Int_t* firstTr=new Int_t[nTrksOrig];
+  Int_t* lastTr=new Int_t[nTrksOrig];
+
+  firstTr[0]=0;
+  for(Int_t iTr=0; iTr<nTrksOrig-1; iTr++){
+    Double_t distz=zTr[iTr+1]-zTr[iTr];
+    Double_t errdistz=TMath::Sqrt(err2zTr[iTr+1]+err2zTr[iTr]);
+    if(errdistz<=0.000001) errdistz=0.000001;
+    if(distz>fDeltaZCutForCluster || (distz/errdistz)>fnSigmaZCutForCluster){
+      lastTr[nClusters]=iTr;
+      firstTr[nClusters+1]=iTr+1;
+      nClusters++;
+    }
+  }
+  lastTr[nClusters]=nTrksOrig-1;
+
+  // Compute vertices
+  AliDebug(1,Form("Number of found clusters %d",nClusters+1));
+  Int_t nFoundVertices=0;
+
+  if (!fMVVertices) fMVVertices = new TObjArray(nClusters+1);
+
+  fMVVertices->Clear();
+  TObjArray cluTrackArr(nTrksOrig);
+  UShort_t *idTrkClu=new UShort_t[nTrksOrig];
+  //  Int_t maxContr=0;
+  //  Int_t maxPos=-1;
+
+  for(Int_t iClu=0; iClu<=nClusters; iClu++){
+    Int_t nCluTracks=lastTr[iClu]-firstTr[iClu]+1;
+    cluTrackArr.Clear();
+    AliDebug(1,Form(" Vertex #%d tracks %d first tr %d  last track %d",iClu,nCluTracks,firstTr[iClu],lastTr[iClu]));
+    Int_t nSelTr=0;
+    for(Int_t iTr=firstTr[iClu]; iTr<=lastTr[iClu]; iTr++){
+      AliExternalTrackParam* t=(AliExternalTrackParam*)trkArrayOrig->At(posOrig[iTr]);
+      if(TMath::Abs(t->GetZ()-zTr[iTr])>0.00001){
+       AliError(Form("Clu %d Track %d zTrack=%f  zVec=%f\n",iClu,iTr,t->GetZ(),zTr[iTr]));
+      }
+      cluTrackArr.AddAt(t,nSelTr);
+      idTrkClu[nSelTr]=idOrig[iTr];
+      AliDebug(1,Form("   Add track %d: id %d, z=%f",iTr,idOrig[iTr],zTr[iTr]));
+      nSelTr++;
+    }
+    AliESDVertex* vert=FindPrimaryVertex(&cluTrackArr,idTrkClu);
+    AliDebug(1,Form("Found vertex in z=%f with %d contributors",vert->GetZv(),
+                vert->GetNContributors()));
+
+    fCurrentVertex=0;
+    if(vert->GetNContributors()>0){
+      nFoundVertices++;
+      fMVVertices->AddLast(vert);
+    }
+    //    if(vert->GetNContributors()>maxContr){
+    //      maxContr=vert->GetNContributors();
+    //      maxPos=nFoundVertices-1;
+    //    }
+  }
+
+  AliDebug(1,Form("Number of found vertices %d (%d)",nFoundVertices,fMVVertices->GetEntriesFast()));
+  // if(maxPos>=0 && maxContr>0){
+  //   AliESDVertex* vtxMax=(AliESDVertex*)fMVVertices->At(maxPos);
+  //   if(fCurrentVertex) delete fCurrentVertex; 
+  //   fCurrentVertex=new AliESDVertex(*vtxMax);
+  // }
+
+  delete [] firstTr;
+  delete [] lastTr;
+  delete [] idTrkClu;
+  delete [] posOrig;
+
+  return;
+
+}