From 30da2335b27961def2ff7ca63e127d8d71339a71 Mon Sep 17 00:00:00 2001 From: masera Date: Thu, 20 Feb 2014 17:27:44 +0100 Subject: [PATCH] New vertexer with IB tracklets --- ITS/UPGRADE/AliITSUClusterLines.cxx | 169 +++++++ ITS/UPGRADE/AliITSUClusterLines.h | 49 ++ ITS/UPGRADE/AliITSUVertexer.cxx | 615 ++++++++++++++++++++++++-- ITS/UPGRADE/AliITSUVertexer.h | 94 +++- ITS/UPGRADE/CMakelibITSUpgradeRec.pkg | 1 + ITS/UPGRADE/ITSUpgradeRecLinkDef.h | 1 + 6 files changed, 886 insertions(+), 43 deletions(-) create mode 100644 ITS/UPGRADE/AliITSUClusterLines.cxx create mode 100644 ITS/UPGRADE/AliITSUClusterLines.h diff --git a/ITS/UPGRADE/AliITSUClusterLines.cxx b/ITS/UPGRADE/AliITSUClusterLines.cxx new file mode 100644 index 00000000000..e6d59004f35 --- /dev/null +++ b/ITS/UPGRADE/AliITSUClusterLines.cxx @@ -0,0 +1,169 @@ +#include "AliITSUClusterLines.h" +#include + +ClassImp(AliITSUClusterLines); + +////////////////////////////////////////////////////////////////////// +// This class is used by the AliITSUVertexer to compute and store // +// information about vertex candidates. // +// The cluster is seeded starting from two lines, then it is // +// possible to attach other lines. Whenever a new line is attached // +// the weight and coefficient matrices are computed. // +// Origin puccio@to.infn.it Feb. 20 2014 // +////////////////////////////////////////////////////////////////////// + + +//____________________________________________________________________________________________________________ +AliITSUClusterLines::AliITSUClusterLines() : TObject(), + fA(), + fB(), + fLabels(), + fV() { + // Default Constructor +} + +//____________________________________________________________________________________________________________ +AliITSUClusterLines::AliITSUClusterLines(UInt_t first, AliStrLine *line1, UInt_t second, AliStrLine *line2,Bool_t weight) : TObject(), + fA(), + fB(), + fLabels(), + fV() { + // Standard constructor + fLabels.push_back(first); + fLabels.push_back(second); + Double_t wmat1[9],wmat2[9],p1[3],p2[3],cd1[3],cd2[3],sq1[3]={1.,1.,1.},sq2[3]={1.,1.,1.}; + line1->GetWMatrix(wmat1); + line2->GetWMatrix(wmat2); + for(Int_t i=0;i<9;++i) fW[i]=wmat1[i]+wmat2[i]; + line1->GetP0(p1); + line2->GetP0(p2); + line1->GetCd(cd1); + line2->GetCd(cd2); + if(weight) { + line1->GetSigma2P0(sq1); + line2->GetSigma2P0(sq2); + } + + Double_t det1=cd1[2]*cd1[2]*sq1[0]*sq1[1]+cd1[1]*cd1[1]*sq1[0]*sq1[2]+cd1[0]*cd1[0]*sq1[1]*sq1[2]; + Double_t det2=cd2[2]*cd2[2]*sq2[0]*sq2[1]+cd2[1]*cd2[1]*sq2[0]*sq2[2]+cd2[0]*cd2[0]*sq2[1]*sq2[2]; + + fA[0]=(cd1[2]*cd1[2]*sq1[1]+cd1[1]*cd1[1]*sq1[2])/det1+(cd2[2]*cd2[2]*sq2[1]+cd2[1]*cd2[1]*sq2[2])/det2; + fA[1]=-cd1[0]*cd1[1]*sq1[2]/det1-cd2[0]*cd2[1]*sq2[2]/det2; + fA[2]=-cd1[0]*cd1[2]*sq1[1]/det1-cd2[0]*cd2[2]*sq2[1]/det2; + fA[3]=(cd1[2]*cd1[2]*sq1[0]+cd1[0]*cd1[0]*sq1[2])/det1+(cd2[2]*cd2[2]*sq2[0]+cd2[0]*cd2[0]*sq2[2])/det2; + fA[4]=-cd1[1]*cd1[2]*sq1[0]/det1-cd2[1]*cd2[2]*sq2[0]/det2; + fA[5]=(cd1[1]*cd1[1]*sq1[0]+cd1[0]*cd1[0]*sq1[1])/det1+(cd2[1]*cd2[1]*sq2[0]+cd2[0]*cd2[0]*sq2[1])/det2; + + fB[0]=(cd1[1]*sq1[2]*(-cd1[1]*p1[0]+cd1[0]*p1[1])+cd1[2]*sq1[1]*(-cd1[2]*p1[0]+cd1[0]*p1[2]))/det1; + fB[0]+=(cd2[1]*sq2[2]*(-cd2[1]*p2[0]+cd2[0]*p2[1])+cd2[2]*sq2[1]*(-cd2[2]*p2[0]+cd2[0]*p2[2]))/det2; + fB[1]=(cd1[0]*sq1[2]*(-cd1[0]*p1[1]+cd1[1]*p1[0])+cd1[2]*sq1[0]*(-cd1[2]*p1[1]+cd1[1]*p1[2]))/det1; + fB[1]+=(cd2[0]*sq2[2]*(-cd2[0]*p2[1]+cd2[1]*p2[0])+cd2[2]*sq2[0]*(-cd2[2]*p2[1]+cd2[1]*p2[2]))/det2; + fB[2]=(cd1[0]*sq1[1]*(-cd1[0]*p1[2]+cd1[2]*p1[0])+cd1[1]*sq1[0]*(-cd1[1]*p1[2]+cd1[2]*p1[1]))/det1; + fB[2]+=(cd2[0]*sq2[1]*(-cd2[0]*p2[2]+cd2[2]*p2[0])+cd2[1]*sq2[0]*(-cd2[1]*p2[2]+cd2[2]*p2[1]))/det2; + + this->ComputeClusterCentroid(); + +} + +//____________________________________________________________________________________________________________ +AliITSUClusterLines::~AliITSUClusterLines() { + // Destructor +} + +//____________________________________________________________________________________________________________ +void AliITSUClusterLines::Add(UInt_t label, AliStrLine *line, Bool_t weight) { + // Add a line to the cluster. It changes the weight matrix of the cluster and its parameters + fLabels.push_back(label); + Double_t wmat[9],p[3],cd[3],sq[3]={1.,1.,1.}; + line->GetWMatrix(wmat); + line->GetP0(p); + line->GetCd(cd); + + for(Int_t i=0;i<9;++i) fW[i]+=wmat[i]; + if(weight) line->GetSigma2P0(sq); + + Double_t det=cd[2]*cd[2]*sq[0]*sq[1]+cd[1]*cd[1]*sq[0]*sq[2]+cd[0]*cd[0]*sq[1]*sq[2]; + fA[0]+=(cd[2]*cd[2]*sq[1]+cd[1]*cd[1]*sq[2])/det; + fA[1]+=-cd[0]*cd[1]*sq[2]/det; + fA[2]+=-cd[0]*cd[2]*sq[1]/det; + fA[3]+=(cd[2]*cd[2]*sq[0]+cd[0]*cd[0]*sq[2])/det; + fA[4]+=-cd[1]*cd[2]*sq[0]/det; + fA[5]+=(cd[1]*cd[1]*sq[0]+cd[0]*cd[0]*sq[1])/det; + + fB[0]+=(cd[1]*sq[2]*(-cd[1]*p[0]+cd[0]*p[1])+cd[2]*sq[1]*(-cd[2]*p[0]+cd[0]*p[2]))/det; + fB[1]+=(cd[0]*sq[2]*(-cd[0]*p[1]+cd[1]*p[0])+cd[2]*sq[0]*(-cd[2]*p[1]+cd[1]*p[2]))/det; + fB[2]+=(cd[0]*sq[1]*(-cd[0]*p[2]+cd[2]*p[0])+cd[1]*sq[0]*(-cd[1]*p[2]+cd[2]*p[1]))/det; + + this->ComputeClusterCentroid(); +} + +//____________________________________________________________________________________________________________ +Int_t AliITSUClusterLines::Compare(const TObject* obj) const { + // Comparison criteria between two clusters + const AliITSUClusterLines *cl=(const AliITSUClusterLines*)obj; + if(fLabels.size()GetSize()) return 1; + if(fLabels.size()>cl->GetSize()) return -1; + return 0; +} + +//____________________________________________________________________________________________________________ +void AliITSUClusterLines::ComputeClusterCentroid() { + // Calculation of the centroid + Double_t *a=fA; + Double_t *b=fB; + Double_t *v=fV; + + Double_t det=a[0]*(a[3]*a[5]-a[4]*a[4])-a[1]*(a[1]*a[5]-a[4]*a[2])+a[2]*(a[1]*a[4]-a[2]*a[3]); + + if(det==0) { + cout << "Could not invert weight matrix" << endl; + return; + } + v[0]=-(b[0]*(a[3]*a[5]-a[4]*a[4])-a[1]*(b[1]*a[5]-a[4]*b[2])+a[2]*(b[1]*a[4]-b[2]*a[3]))/det; + v[1]=-(a[0]*(b[1]*a[5]-b[2]*a[4])-b[0]*(a[1]*a[5]-a[4]*a[2])+a[2]*(a[1]*b[2]-a[2]*b[1]))/det; + v[2]=-(a[0]*(a[3]*b[2]-b[1]*a[4])-a[1]*(a[1]*b[2]-b[1]*a[2])+b[0]*(a[1]*a[4]-a[2]*a[3]))/det; +} + +//____________________________________________________________________________________________________________ +void AliITSUClusterLines::GetCovMatrix(Float_t cov[6]) { + // Returns the covariance matrix (single precision) + Double_t *w=fW; + Double_t den=w[0]*(w[4]*w[8]-w[5]*w[7])-w[1]*(w[3]*w[8]-w[5]*w[6])+w[2]*(w[3]*w[7]-w[4]*w[6]); + if(den==0) { + cout << "Could not invert weight matrix" << endl; + return; + } + cov[0]=(w[4]*w[8]-w[5]*w[7])/den; + cov[1]=-(w[1]*w[8]-w[7]*w[2])/den; + cov[2]=(w[1]*w[5]-w[4]*w[2])/den; + cov[3]=(w[0]*w[8]-w[6]*w[2])/den; + cov[4]=-(w[0]*w[5]-w[3]*w[2])/den; + cov[5]=(w[0]*w[4]-w[1]*w[3])/den; +} + +//____________________________________________________________________________________________________________ +void AliITSUClusterLines::GetCovMatrix(Double_t cov[6]) { + // Returns the covariance matrix (double precision) + Double_t *w=fW; + Double_t den=w[0]*(w[4]*w[8]-w[5]*w[7])-w[1]*(w[3]*w[8]-w[5]*w[6])+w[2]*(w[3]*w[7]-w[4]*w[6]); + if(den==0) { + cout << "Could not invert weight matrix" << endl; + return; + } + cov[0]=(w[4]*w[8]-w[5]*w[7])/den; + cov[1]=-(w[1]*w[8]-w[7]*w[2])/den; + cov[2]=(w[1]*w[5]-w[4]*w[2])/den; + cov[3]=(w[0]*w[8]-w[6]*w[2])/den; + cov[4]=-(w[0]*w[5]-w[3]*w[2])/den; + cov[5]=(w[0]*w[4]-w[1]*w[3])/den; + +} + +//____________________________________________________________________________________________________________ +Bool_t AliITSUClusterLines::IsEqual(const TObject* obj) const { + // Comparison criteria between two clusters + const AliITSUClusterLines *cl=(const AliITSUClusterLines*)obj; + if(fLabels.size()==cl->GetSize()) return kTRUE; + return kFALSE; +} + diff --git a/ITS/UPGRADE/AliITSUClusterLines.h b/ITS/UPGRADE/AliITSUClusterLines.h new file mode 100644 index 00000000000..50586746ff9 --- /dev/null +++ b/ITS/UPGRADE/AliITSUClusterLines.h @@ -0,0 +1,49 @@ +#ifndef ALIITSUCLUSTERLINES_H +#define ALIITSUCLUSTERLINES_H +/* Copyright(c) 2009-2014, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +///////////////////////////////////////////////////////////////////////////// +// Class intended to gather and compute the parameters of vertex candidate // +///////////////////////////////////////////////////////////////////////////// + +#include +#include +#include "AliStrLine.h" + + +class AliITSUClusterLines : public TObject { + public: + AliITSUClusterLines(); + AliITSUClusterLines(UInt_t first, AliStrLine *firstL, UInt_t second, AliStrLine *secondL,Bool_t=kFALSE); + virtual ~AliITSUClusterLines(); + + void Add(UInt_t label, AliStrLine *line, Bool_t weight=kFALSE); + void ComputeClusterCentroid(); + inline UInt_t GetSize() const { return fLabels.size(); } + inline Int_t* GetLabels(UInt_t &n) {n=fLabels.size(); return &fLabels[0]; } + inline void GetA(Float_t a[3]) { for(Int_t i=0; i<3; ++i) a[i]=fA[i]; } + inline void GetA(Double_t a[3]) { for(Int_t i=0; i<3; ++i) a[i]=fA[i]; } + inline void GetB(Float_t b[3]) { for(Int_t i=0; i<3; ++i) b[i]=fB[i]; } + inline void GetB(Double_t b[3]) { for(Int_t i=0; i<3; ++i) b[i]=fB[i]; } + void GetCovMatrix(Float_t cov[6]); + void GetCovMatrix(Double_t cov[6]); + inline void GetVertex(Float_t p[3]) { for(Int_t i=0; i<3; ++i) p[i]=fV[i]; } + inline void GetVertex(Double_t p[3]) { for(Int_t i=0; i<3; ++i) p[i]=fV[i]; } + inline void GetWeight(Float_t w[9]) { for(Int_t i=0; i<9; ++i) w[i]=fW[i]; } + inline void GetWeight(Double_t w[9]) { for(Int_t i=0; i<9; ++i) w[i]=fW[i]; } + // + virtual Bool_t IsSortable() const {return kTRUE;} + virtual Bool_t IsEqual(const TObject* obj) const; + virtual Int_t Compare(const TObject* obj) const; + + protected: + Double_t fA[6]; // AX=B + Double_t fB[3]; // AX=B + vector fLabels; // labels + Double_t fV[3]; // vertex candidate + Double_t fW[9]; // weight matrix + ClassDef(AliITSUClusterLines,1); +}; + +#endif diff --git a/ITS/UPGRADE/AliITSUVertexer.cxx b/ITS/UPGRADE/AliITSUVertexer.cxx index 7d0604f540e..c79c6137850 100644 --- a/ITS/UPGRADE/AliITSUVertexer.cxx +++ b/ITS/UPGRADE/AliITSUVertexer.cxx @@ -1,40 +1,163 @@ -/************************************************************************** - * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. * - * * - * Author: The ALICE Off-line Project. * - * Contributors are mentioned in the code where appropriate. * - * * - * Permission to use, copy, modify and distribute this software and its * - * documentation strictly for non-commercial purposes is hereby granted * - * without fee, provided that the above copyright notice appears in all * - * copies and that both the copyright notice and this permission notice * - * appear in the supporting documentation. The authors make no claims * - * about the suitability of this software for any purpose. It is * - * provided "as is" without express or implied warranty. * - **************************************************************************/ -#include "AliITSUVertexer.h" +#include +#include +#include +#include +#include #include "AliESDVertex.h" +#include "AliITSUClusterPix.h" +#include "AliITSUVertexer.h" +#include "AliStrLine.h" +#include "AliVertexerTracks.h" +#include "AliITSUClusterLines.h" #include "AliLog.h" -#include -#include - -///////////////////////////////////////////////////////////////// -// THIS IS A DUMMY CLASS -//////////////////////////////////////////////////////////////// -ClassImp(AliITSUVertexer) +using TMath::Abs; +using TMath::Sqrt; +using TMath::ATan2; +using TMath::TwoPi; +using TMath::BubbleLow; +////////////////////////////////////////////////////////////////////// +// This class is used to compute the position of all the primary // +// vertices in a single event using the upgraded ITS. // +// Optimizations ongoing. // +// Origin puccio@to.infn.it Feb. 20 2014 // +////////////////////////////////////////////////////////////////////// - -//______________________________________________________________________ -AliITSUVertexer::AliITSUVertexer():AliVertexer() +//_____________________________________________________________________________________________ +AliITSUVertexer::AliITSUVertexer(Double_t phicut, Double_t zcut, Double_t paircut,Double_t clcut, Int_t cclcut) : AliVertexer(), + fClusterContribCut(cclcut), + fClusterCut(clcut), + fClusterIndex(), + fClusterPhi(), + fClusters(), + fLines("AliStrLine",1000), + fLinesClusters("AliITSUClusterLines.h",1000), + fLinesPhi(0), + fNoClusters(0), + fNoLines(0), + fNoVertices(0), + fPairCut(paircut), + fPhiCut(phicut), + fZCut(zcut), + fUsedClusters(), + fUsedLines(), + fVertices() +#ifdef MC_CHECK +,fGoodLines(0),fGoodLinesPhi(0),fParticleId(0) +#endif { + // Standard I/O constructor } - -//______________________________________________________________________ +//_____________________________________________________________________________________________ AliITSUVertexer::~AliITSUVertexer() { // Destructor + Reset(); +} + +//_____________________________________________________________________________________________ +void AliITSUVertexer::FindVerticesForCurrentEvent() { + // Try to find all the primary vertices in the current + fNoVertices=0; + FindTracklets(); + if(fNoLines<2) { + fVertices.push_back(AliESDVertex()); + return;// AliESDVertex(); + } + + // fVertices.push_back(AliVertexerTracks::TrackletVertexFinder(&fLines,1)); + //fNoVertices=1; + fUsedLines=new Short_t[fNoLines]; + for(UInt_t i=0;iGetDCA(line2)<=fPairCut) { + //cout << fNoClusters <<" " << i1 << " " << i2 << " "; + new(fLinesClusters[fNoClusters])AliITSUClusterLines(i1,line1,i2,line2); + AliITSUClusterLines* current=(AliITSUClusterLines*)fLinesClusters.At(fNoClusters); + Double_t p[3]; + current->GetVertex(p); + if((p[0]*p[0]+p[1]*p[1])>=4) { // Beam pipe check + fLinesClusters.RemoveAt(fNoClusters); + fLinesClusters.Compress(); + break; + } + fUsedLines[i1]=fNoClusters; + fUsedLines[i2]=fNoClusters; + for(UInt_t i3=0;i3PrintStatus(); + if(line3->GetDistFromPoint(p)<=fPairCut) { + //cout << i3 << " "; + current->Add(i3,line3); + fUsedLines[i3]=fNoClusters; + current->GetVertex(p); + } + //cout << "Ballo la ula" << endl; + } + ++fNoClusters; + //cout << endl; + break; + } + } + } + + fLinesClusters.Sort(); + + for(UInt_t i0=0;i0GetVertex(p0); + for(UInt_t i1=i0+1;i1GetVertex(p1); + if (TMath::Abs(p0[2]-p1[2])<=fClusterCut) { + Double_t distance=(p0[0]-p1[0])*(p0[0]-p1[0])+(p0[1]-p1[1])*(p0[1]-p1[1])+(p0[2]-p1[2])*(p0[2]-p1[2]); + //Bool_t flag=kFALSE; + if(distance<=fPairCut*fPairCut) { + UInt_t n=0; + Int_t *labels=clu1->GetLabels(n); + for(UInt_t icl=0; iclAdd(labels[icl],(AliStrLine*)fLines.At(labels[icl])); + clu0->GetVertex(p0); + //flag=kTRUE; + } + fLinesClusters.RemoveAt(i1); + fLinesClusters.Compress(); + fNoClusters--; + i1--; + //if(flag) i1=10; + } + } + } + + for(Int_t i0=fNoClusters-1;i0>=0; --i0) { + AliITSUClusterLines *clu0 = (AliITSUClusterLines*)fLinesClusters.At(i0); + Int_t size=clu0->GetSize(); + if(size1) { + fLinesClusters.RemoveAt(i0); + fLinesClusters.Compress(); + fNoClusters--; + continue; + } + Double_t p0[3],cov[6]; + clu0->GetVertex(p0); + clu0->GetCovMatrix(cov); + if((p0[0]*p0[0]+p0[1]*p0[1])<1.98*1.98) { + fVertices.push_back(AliESDVertex(p0,cov,99999.,size)); + } + } + fNoVertices=fVertices.size(); + + return;// AliVertexerTracks::TrackletVertexFinder(&fLines,0); } //______________________________________________________________________ @@ -49,3 +172,441 @@ AliESDVertex* AliITSUVertexer::FindVertexForCurrentEvent(TTree *) return vtx; } +//_____________________________________________________________________________________________ +Int_t AliITSUVertexer::MatchPoints(UShort_t layer, Double_t anchor, Double_t *p0, Double_t *p1) { + // Method for matching clusters with similar phi and z inside the range [zmin,zmax] + AliDebug(3,Form("Matching points on layer %i...\n",layer)); + + Double_t xyz[3][3]; // {{0: 'min', 1: 'mid', 2: 'max'},{0: 'x', 1: 'y', 2: 'z'}} + AliITSUClusterPix* cl[3]; // {0: 'min', 1: 'mid', 2: 'max'} + Double_t phi[3]; // {0: 'min', 1: 'mid', 2: 'max'} + + Int_t nocl=fClusters[layer]->GetEntriesFast(); + Int_t imin=0,imax=nocl!=0 ? nocl-1 : nocl,imid=(imax+imin)/2; + + Double_t a=0,b=0,c=0; + if(layer==2) { + a=(p0[0]-p1[0])*(p0[0]-p1[0])+(p0[1]-p1[1])*(p0[1]-p1[1]); + b=(p1[0]-p0[0])*p0[0]+(p1[1]-p0[1])*p0[1]; + c=p0[0]*p0[0]+p0[1]*p0[1]; + } + + Int_t flag=-1; + + while (imax > imin) { + while(fUsedClusters[layer][fClusterIndex[layer][imin]]==kTRUE&&imin0) --imax; + if(imaxAt(fClusterIndex[layer][imin]); + cl[0]->GetGlobalXYZ(xyz[0]); + phi[0] = fClusterPhi[layer][fClusterIndex[layer][imin]]; + AliDebug(4,Form("Cluster min: %i, %i, %f, %i\n",imin,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][imin]])); + if((Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)) { + AliDebug(4,Form("Ok min: %i \n",imin)); + if(layer==2) { + if(fUsedClusters[layer][fClusterIndex[layer][imin]]) { + flag=1; + break; + } + c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]); + Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a; + Double_t zmin=z-fZCut,zmax=z+fZCut; + if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) { + AliDebug(4,Form("Ok Z: %i \n",imin)); + return fClusterIndex[layer][imin]; + } + AliDebug(4,Form("No Z match: %i \n",imin)); + c=p0[0]*p0[0]+p0[1]*p0[1]; + flag=1; + break; + } else if(!fUsedClusters[layer][fClusterIndex[layer][imin]]) return fClusterIndex[layer][imin]; + } + + cl[2] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imax]); + cl[2]->GetGlobalXYZ(xyz[2]); + phi[2] = fClusterPhi[layer][fClusterIndex[layer][imax]]; + AliDebug(4,Form("Cluster max: %i, %i, %f, %i\n",imax,cl[2]->GetLabel(0),phi[2]-anchor,fUsedClusters[layer][fClusterIndex[layer][imax]])); + if((Abs(phi[2]-anchor)<=fPhiCut||Abs(Abs(phi[2]-anchor)-TMath::TwoPi())<=fPhiCut)) { + AliDebug(4,Form("Ok max: %i \n",imax)); + if(layer==2) { + if(fUsedClusters[layer][fClusterIndex[layer][imax]]) { + flag=3; + break; + } + c-=(xyz[2][0]*xyz[2][0]+xyz[2][1]*xyz[2][1]); + Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a; + Double_t zmin=z-fZCut,zmax=z+fZCut; + if(zmin<=xyz[2][2]&&xyz[2][2]<=zmax) { + AliDebug(4,Form("Ok Z: %i \n",imax)); + return fClusterIndex[layer][imax]; + } + c=p0[0]*p0[0]+p0[1]*p0[1]; + AliDebug(4,Form("No Z match: %i \n",imax)); + flag=3; + break; + } else if(!fUsedClusters[layer][fClusterIndex[layer][imax]]) return fClusterIndex[layer][imax]; + } + + imid=(imax+imin)/2; + if(imid==imin) return -1; + Int_t step=1,sign=-1; + while(fUsedClusters[layer][fClusterIndex[layer][imid]]) { + imid=imid+step; + sign*=-1; + step=(step+sign)*(-1); + if(imid>=imax||imid<=imin) return -1; + } + + cl[1] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imid]); + cl[1]->GetGlobalXYZ(xyz[1]); + phi[1] = fClusterPhi[layer][fClusterIndex[layer][imid]]; + AliDebug(4,Form("Cluster mid: %i, %i, %f, %i \n",imid,cl[1]->GetLabel(0),phi[1]-anchor,fUsedClusters[layer][fClusterIndex[layer][imid]])); + //cout << imin << " " << imid << " " << imax << endl; + //cout << fClusterIndex[layer][imid] << endl; + if((Abs(phi[1]-anchor)<=fPhiCut||Abs(Abs(phi[1]-anchor)-TMath::TwoPi())<=fPhiCut)) { + AliDebug(4,Form("Ok mid: %i \n",imid)); + if(layer==2) { + c-=(xyz[1][0]*xyz[1][0]+xyz[1][1]*xyz[1][1]); + Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a; + Double_t zmin=z-fZCut,zmax=z+fZCut; + if(zmin<=xyz[1][2]&&xyz[1][2]<=zmax) { + AliDebug(4,Form("Ok Z: %i \n",imid)); + return fClusterIndex[layer][imid]; + } + AliDebug(4,Form("No Z match: %i \n",imid)); + c=p0[0]*p0[0]+p0[1]*p0[1]; + flag=2; + break; + } else if(!fUsedClusters[layer][fClusterIndex[layer][imid]]) return fClusterIndex[layer][imid]; + } + + // determine which subarray to search + if (phi[1] < anchor) { + // change min index to search upper subarray + AliDebug(4,"Case minor\n"); + imin = imid + 1; + --imax; + } else if (phi[1] > anchor) { + AliDebug(4,"Case greater\n"); + // change max index to search lower subarray + ++imin; + imax = imid - 1; + } else if(!fUsedClusters[layer][fClusterIndex[layer][imid]]) { + return fClusterIndex[layer][imid]; + } else return -1; + } + + if(flag>-1) { + AliDebug(4,"Flag issued, starting forward backward check\n"); + Int_t start=imid; + switch(flag) { + case 1: + start=imin; + break; + case 2: + start=imid; + break; + case 3: + start=imax; + break; + } + + Int_t curr=start-1; + Bool_t lap=kFALSE; + while(1){ + if(curr==-1&&!lap) { + lap=kTRUE; + curr=nocl-1; + } else if(lap) break; + cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][curr]); + cl[0]->GetGlobalXYZ(xyz[0]); + phi[0] = fClusterPhi[layer][fClusterIndex[layer][curr]]; + AliDebug(4,Form("Looking backward: %i, %i, %f, %i\n",curr,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][curr]])); + if(!(Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)) break; + if(fUsedClusters[layer][fClusterIndex[layer][curr]]) { + curr--; + continue; + } + + AliDebug(4,Form("Ok backward: %i \n",curr)); + + c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]); + Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a; + Double_t zmin=z-fZCut,zmax=z+fZCut; + if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) { + AliDebug(4,Form("Ok Z: %i \n",curr)); + return fClusterIndex[layer][curr]; + } + AliDebug(4,Form("No Z match: %i \n",curr)); + c=p0[0]*p0[0]+p0[1]*p0[1]; + curr--; + } + + lap=kFALSE; + curr=start+1; + while(1){ + if(curr==nocl&&!lap) { + curr=0; + lap=kTRUE; + } else if(lap) break; + cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][curr]); + cl[0]->GetGlobalXYZ(xyz[0]); + phi[0] = fClusterPhi[layer][fClusterIndex[layer][curr]]; + AliDebug(4,Form("Looking forward: %i, %i, %f, %i\n",curr,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][curr]])); + if(!(Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)) break; + if(fUsedClusters[layer][fClusterIndex[layer][curr]]) { + curr++; + continue; + } + AliDebug(4,Form("Ok forward: %i \n",curr)); + c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]); + Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a; + Double_t zmin=z-fZCut,zmax=z+fZCut; + if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) { + AliDebug(4,Form("Ok Z: %i \n",curr)); + return fClusterIndex[layer][curr]; + } + AliDebug(4,Form("No Z match: %i \n",curr)); + c=p0[0]*p0[0]+p0[1]*p0[1]; + curr++; + } + } + + if(imax==imin&&imax!=0) { + cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imin]); + cl[0]->GetGlobalXYZ(xyz[0]); + phi[0] = fClusterPhi[layer][fClusterIndex[layer][imin]]; + AliDebug(4,Form("Cluster eq: %i, %i, %f, %i\n",imin,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][imin]])); + if((Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)&&!fUsedClusters[layer][fClusterIndex[layer][imin]]) { + AliDebug(4,Form("Ok eq: %i \n",imin)); + if(layer==2) { + c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]); + Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a; + Double_t zmin=z-fZCut,zmax=z+fZCut; + if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) { + AliDebug(4,Form("Ok Z eq: %i \n",imin)); + return fClusterIndex[layer][imin]; + } + AliDebug(4,Form("No Z eq: %i \n",imin)); + c=p0[0]*p0[0]+p0[1]*p0[1]; + } else return fClusterIndex[layer][imin]; + } + } + // no match found :( + return -1; + +} + +//_____________________________________________________________________________________________ +void AliITSUVertexer::PrintStatus() const { + // Prints all the cuts and important data members for the current status + cout << "Cut on phi: " << fPhiCut << endl; + cout << "Cut on z: " << fZCut << endl; +} + +//_____________________________________________________________________________________________ +void AliITSUVertexer::Reset() { + // Resets the vertexer for a new event (or for its destruction) + AliDebug(2,"Resetting the vertexer...\n"); + fNoVertices=0; + for(Int_t i=0;i<3;++i) { + delete[] fUsedClusters[i]; + delete[] fClusterIndex[i]; + delete[] fClusterPhi[i]; + } + if(fNoLines>2) { + delete []fUsedLines; + } + + fLinesPhi=0; + fLines.Clear(); + fLinesClusters.Clear(); + fVertices.clear(); + + #ifdef MC_CHECK + fGoodLines=0; + fGoodLinesPhi=0; + delete[] fParticleId; + #endif + + //delete fVertices; +} + +//_____________________________________________________________________________________________ +void AliITSUVertexer::FindTracklets() { + // It combines recpoints over the first three layers to define a list of tracklets + // Here it uses the ordered list of points. Second and third layers are ordered using phi + AliDebug(2,"Calling the trackleter...\n"); + + UInt_t noPntL[3]; + for(UShort_t i=0;i<3;++i) { + noPntL[i]=fClusters[i]->GetEntries(); + fUsedClusters[i]=new Bool_t[noPntL[i]]; + for(UInt_t ii=0;iiAt(i0); + cluster0->GetGlobalXYZ(p0); + vector tmp; + Int_t i1=0; + Int_t label0=cluster0->GetLabel(0); + AliDebug(4,Form("Layer 0: %i\n",label0)); + while(i1>=0) { + i1 = MatchPoints(1,fClusterPhi[0][i0]); + if(i1<0) break; + + ++nolinesphi; + AliITSUClusterPix* cluster1=(AliITSUClusterPix*)fClusters[1]->At(i1); + cluster1->GetGlobalXYZ(p1); + tmp.push_back(i1); + fUsedClusters[1][i1]=kTRUE; + + Int_t i2 = MatchPoints(2,fClusterPhi[1][i1],p0,p1); + if(i2<0) continue; + + #ifdef MC_CHECK + CheckMC(i0,i1,i2); + #endif + + AliITSUClusterPix* cluster2=(AliITSUClusterPix*)fClusters[2]->At(i2); + cluster2->GetGlobalXYZ(pp2); + + fUsedClusters[0][i0]=kTRUE; + fUsedClusters[2][i2]=kTRUE; + + // Errors to be checked... + + Float_t cov0[6],cov1[6],cov2[6]; + cluster0->GetGlobalCov(cov0); + cluster1->GetGlobalCov(cov1); + cluster2->GetGlobalCov(cov2); + + //Error on tracklet direction near the vertex + Double_t rad1=TMath::Sqrt(p0[0]*p0[0]+p0[1]*p0[1]); + Double_t rad2=TMath::Sqrt(p1[0]*p1[0]+p1[1]*p1[1]); + Double_t factor=(rad1+rad2)/(rad2-rad1); + + //Curvature error + Double_t curvErr=0; + Double_t bField=0.5; + Double_t meanPtSelTrk=0.630; + Double_t curvRadius=meanPtSelTrk/(0.3*bField)*100; //cm + Double_t dRad=TMath::Sqrt((p0[0]-p1[0])*(p0[0]-p1[0])+(p0[1]-p1[1])*(p0[1]-p1[1])); + Double_t aux=dRad/2.+rad1; + curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm + + Double_t sq[3],wmat[9]={1,0,0,0,1,0,0,0,1}; + sq[0]=(cov0[0]+curvErr*curvErr/2.)*factor*factor;//cov1[0]+cov2[0]); + sq[1]=(cov0[3]+curvErr*curvErr/2.)*factor*factor;//cov1[3]+cov2[3]); + sq[2]=(cov0[5])*factor*factor;//cov1[5]+cov2[5]); + + // Multiple scattering + Double_t meanPSelTrk=0.875; + Double_t pOverMass=meanPSelTrk/0.140; + Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass); + Double_t p2=meanPSelTrk*meanPSelTrk; + Double_t rBP=1.98; // Beam pipe radius + Double_t dBP=0.08/35.3; // 800 um of Be + Double_t dL1=0.01; //approx. 1% of radiation length + Double_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP; + Double_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1; + Double_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1)); + Double_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP)); + for(Int_t ico=0; ico<3;ico++){ + sq[ico]+=rtantheta1*rtantheta1*factor*factor/3.; + sq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.; + } + + if(sq[0]!=0) wmat[0]=1/sq[0]; + if(sq[1]!=0) wmat[4]=1/sq[1]; + if(sq[2]!=0) wmat[8]=1/sq[2]; + + /*Int_t label0=cluster0->GetLabel(0); + Int_t label1=cluster1->GetLabel(0); + Int_t label2=cluster2->GetLabel(0);*/ + //cout << label0 << " " << label1 << " "<Exec(Form("echo %i >> trkl",cluster0->GetLabel(0)+cluster1->GetLabel(0))); + + new(fLines[fNoLines++])AliStrLine(p0,sq,wmat,p1,kTRUE); + break; + } + for(UInt_t itmp=0; itmpPrintStatus(); + //printf("(%f,%f,%f) and (%f,%f,%f)\n",p0[0],p0[1],p0[2],p1[0],p1[1],p1[2]); + } + + fLinesPhi=nolinesphi; + //cout << Form("=======================================================\nNolinesphi: %i, nolines: %i\nGood Nolinesphi: %i, good nolines: %i\n",nolinesphi,fNoLines,fGoodLinesPhi,fGoodLines); +} + +//___________________________________________________________________________ +void AliITSUVertexer::SetClusters(TClonesArray *clr, const UShort_t i) { + // Reading of the clusters on the first three layer of the upgraded ITS + fClusters[i]=clr; + Int_t nocl=clr->GetEntriesFast(); + if(nocl==0) { + fClusterPhi[i]=new Double_t[1]; + fClusterPhi[i][0]=-999999; + fClusterIndex[i]=new Int_t[1]; + fClusterIndex[i][0]=0; + } else { + fClusterPhi[i]=new Double_t[nocl]; + fClusterIndex[i]=new Int_t[nocl]; + for(Int_t k=0;kAt(k); + Double_t pt[3]; + cl->GetGlobalXYZ(pt); + fClusterPhi[i][k]=ATan2(pt[1],pt[0]); + } + BubbleLow(nocl,fClusterPhi[i],fClusterIndex[i]); + } +} + +#ifdef MC_CHECK +//___________________________________________________________________________ +Bool_t AliITSUVertexer::CheckMC(UInt_t i0, UInt_t i1, UInt_t i2) { + // Debugging function + //cout << "Checking MC truth" << endl; + Int_t label=0; + Bool_t flag=kFALSE; + AliITSUClusterPix* p0=(AliITSUClusterPix*)fClusters[0]->At(i0); + AliITSUClusterPix* p1=(AliITSUClusterPix*)fClusters[1]->At(i1); + for(Int_t i=0; i<3; ++i) { + label=p0->GetLabel(i); + for(Int_t j=0; j<3; ++j) { + if(label==p1->GetLabel(j)&&label>0) { + fGoodLinesPhi++; + flag=kTRUE; + break; + } + } + if(flag) break; + } + + if(!flag) return kFALSE; + + AliITSUClusterPix* p2=(AliITSUClusterPix*)fClusters[2]->At(i2); + for(Int_t j=0; j<3; ++j) { + if(label==p2->GetLabel(j)) { + fParticleId[fGoodLines]=label; + fGoodLines++; + //cout << label << endl; + return kTRUE; + } + } + return kFALSE; +} +#endif diff --git a/ITS/UPGRADE/AliITSUVertexer.h b/ITS/UPGRADE/AliITSUVertexer.h index 1b576e5f80d..8c070ff6a8d 100644 --- a/ITS/UPGRADE/AliITSUVertexer.h +++ b/ITS/UPGRADE/AliITSUVertexer.h @@ -1,31 +1,93 @@ -#ifndef ALIITSUVERTEXERZ_H -#define ALIITSUVERTEXERZ_H +#ifndef ALIITSUVERTEXER_H +#define ALIITSUVERTEXER_H +//#define MC_CHECK // comment out to enable MC checks for debugging -#include +#include "AliVertexer.h" -/////////////////////////////////////////////////////////////////// -// // -// Dummy vertexer plugin // -// // -/////////////////////////////////////////////////////////////////// +/* Copyright(c) 2009-2014, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +///////////////////////////////////////////////////////////////////////////// +// Class for the reconstruction of the primary vertices using ITSU // +///////////////////////////////////////////////////////////////////////////// +class TClonesArray; class AliESDVertex; class AliITSUVertexer : public AliVertexer { - public: - AliITSUVertexer(); + + // Constructors and destructors + AliITSUVertexer(Double_t phicut=0.02,Double_t zcut=0.005,Double_t paircut=0.1, Double_t clustercut=0.5, Int_t clcontrib=3); virtual ~AliITSUVertexer(); - virtual AliESDVertex* FindVertexForCurrentEvent(TTree *itsClusterTree); - void PrintStatus() const {} + + // Public methods + virtual AliESDVertex* GetAllVertices(Int_t& nVert) const { nVert=fNoVertices; return (AliESDVertex*)&fVertices[0]; }; + virtual AliESDVertex* FindVertexForCurrentEvent(TTree *); + void FindVerticesForCurrentEvent(); + virtual void PrintStatus() const; + void Reset(); + + // Getters + UInt_t GetNoLines() const { return fNoLines; } + UShort_t GetNumOfVertices() const { return fNoVertices; } + + + // Setters + void SetClusters(TClonesArray *clr, const UShort_t i); + void SetPhiCut(const Double_t phicut) { fPhiCut=phicut; } + void SetZCut(const Double_t zcut) { fZCut=zcut; } + + #ifdef MC_CHECK + // Debug + MC truth + UInt_t* GetParticleId(UInt_t &num) const { num=fGoodLines; return fParticleId; } + UInt_t GetGoodLines() const { return fGoodLines; } + UInt_t GetGoodLinesPhi() const { return fGoodLinesPhi; } + UInt_t GetLinesPhi() const { return fLinesPhi; } + #endif protected: + // Methods + AliITSUVertexer(AliITSUVertexer&); + AliITSUVertexer& operator=(const AliITSUVertexer& other); + void AddToCluster(UInt_t line,Bool_t weight=kFALSE,Int_t cl=-1); + void CleanAndOrderClusters(); + void Clusterize(UInt_t l1, UInt_t l2, Bool_t weight=kFALSE); + void ComputeClusterCentroid(UInt_t cl); + void FindTracklets(); + Int_t MatchPoints(UShort_t layer, Double_t anchor, Double_t *p0=0x0, Double_t *p1=0x0); + void MoveLabels(Short_t start, Short_t end); + + // Data members + Int_t fClusterContribCut; + Double_t fClusterCut; + Int_t *fClusterIndex[3]; // AliITSUClusterPix index + Double_t *fClusterPhi[3]; // Phi of clusters + TClonesArray *fClusters[3]; //! array of pointers to TClonesArray of AliITSUClusterPix not owned by this class + TClonesArray fLines; //! array of tracklets + TClonesArray fLinesClusters; // array of vertex candidates + UInt_t fLinesPhi; // number of tracklets built by using the first two layers + UInt_t fNoClusters; // number of clusters + UInt_t fNoLines; // number of tracklets + UShort_t fNoVertices; // number of vertices + Double_t fPairCut; // cut on pair + Double_t fPhiCut; // cut on deltaphi for cluster matching among first two layers + Double_t fZCut; // cut on deltatheta for cluster matching among first two layers and the third one + Bool_t *fUsedClusters[3]; // flag for used clusters in tracklet formation + Short_t *fUsedLines; // flag for used lines + vector fVertices; // array of vertices + + #ifdef MC_CHECK + // MC truth methods + Bool_t CheckMC(UInt_t,UInt_t,UInt_t); - private: - AliITSUVertexer(const AliITSUVertexer& vtxr); - AliITSUVertexer& operator=(const AliITSUVertexer& vtxr ); + // MC truth data members + UInt_t fGoodLines; + UInt_t fGoodLinesPhi; + UInt_t *fParticleId; + #endif - ClassDef(AliITSUVertexer,1 ); + ClassDef(AliITSUVertexer,1) }; #endif diff --git a/ITS/UPGRADE/CMakelibITSUpgradeRec.pkg b/ITS/UPGRADE/CMakelibITSUpgradeRec.pkg index 7fbdc440109..afdc6339487 100644 --- a/ITS/UPGRADE/CMakelibITSUpgradeRec.pkg +++ b/ITS/UPGRADE/CMakelibITSUpgradeRec.pkg @@ -30,6 +30,7 @@ set ( SRCS AliITSURecoParam.cxx AliITSUReconstructor.cxx AliITSUClusterizer.cxx + AliITSUClusterLines.cxx AliITSURecoSens.cxx AliITSURecoLayer.cxx AliITSURecoDet.cxx diff --git a/ITS/UPGRADE/ITSUpgradeRecLinkDef.h b/ITS/UPGRADE/ITSUpgradeRecLinkDef.h index e54177b6634..917082fffe3 100644 --- a/ITS/UPGRADE/ITSUpgradeRecLinkDef.h +++ b/ITS/UPGRADE/ITSUpgradeRecLinkDef.h @@ -16,6 +16,7 @@ #pragma link C++ class AliITSURecoParam+; #pragma link C++ class AliITSUReconstructor+; #pragma link C++ class AliITSUClusterizer+; +#pragma link C++ class AliITSUClusterLines+; #pragma link C++ class AliITSURecoSens+; #pragma link C++ class AliITSURecoLayer+; #pragma link C++ class AliITSURecoDet+; -- 2.43.0