X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSVertexer3D.cxx;h=e49a9b7ca80aec8a759b768864620737848a19a4;hb=bc876f8c657e38293f75b7383e583f8246926eb2;hp=421f55f3d67eb1eebd196c460b14a602abfe20fc;hpb=5951029fa093300f404f85e5dbbc6b6583326b99;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSVertexer3D.cxx b/ITS/AliITSVertexer3D.cxx index 421f55f3d67..e49a9b7ca80 100644 --- a/ITS/AliITSVertexer3D.cxx +++ b/ITS/AliITSVertexer3D.cxx @@ -13,12 +13,14 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ #include +#include "AliRunLoader.h" #include "AliESDVertex.h" #include "AliLog.h" #include "AliStrLine.h" #include "AliTracker.h" #include "AliITSDetTypeRec.h" #include "AliITSRecPoint.h" +#include "AliITSRecPointContainer.h" #include "AliITSgeomTGeo.h" #include "AliVertexerTracks.h" #include "AliITSVertexer3D.h" @@ -27,10 +29,12 @@ ///////////////////////////////////////////////////////////////// // this class implements a method to determine // the 3 coordinates of the primary vertex -// for p-p collisions -// It can be used successfully with Pb-Pb collisions +// optimized for +// p-p collisions //////////////////////////////////////////////////////////////// +const Int_t AliITSVertexer3D::fgkMaxNumOfClDefault = 1000; + ClassImp(AliITSVertexer3D) /* $Id$ */ @@ -45,6 +49,7 @@ AliITSVertexer3D::AliITSVertexer3D(): fCutOnPairs(0.), fCoarseMaxRCut(0.), fMaxRCut(0.), + fMaxRCut2(0.), fZCutDiamond(0.), fMaxZCut(0.), fDCAcut(0.), @@ -54,9 +59,13 @@ AliITSVertexer3D::AliITSVertexer3D(): fUsedCluster(kMaxCluPerMod*kNSPDMod), fZHisto(0), fDCAforPileup(0.), + fDiffPhiforPileup(0.), fBinSizeR(0.), fBinSizeZ(0.), - fPileupAlgo(0) + fPileupAlgo(0), + fMaxNumOfCl(fgkMaxNumOfClDefault), + fDoDownScale(kFALSE), + fGenerForDownScale(0) { // Default constructor SetCoarseDiffPhiCut(); @@ -64,6 +73,7 @@ AliITSVertexer3D::AliITSVertexer3D(): SetCutOnPairs(); SetCoarseMaxRCut(); SetMaxRCut(); + SetMaxRCutAlgo2(); SetZCutDiamond(); SetMaxZCut(); SetDCACut(); @@ -71,12 +81,14 @@ AliITSVertexer3D::AliITSVertexer3D(): SetMeanPSelTracks(); SetMeanPtSelTracks(); SetMinDCAforPileup(); + SetDeltaPhiforPileup(); SetPileupAlgo(); SetBinSizeR(); SetBinSizeZ(); - Float_t binsize=0.02; // default 200 micron + Double_t binsize=0.02; // default 200 micron Int_t nbins=static_cast(1+2*fZCutDiamond/binsize); fZHisto=new TH1F("hz","",nbins,-fZCutDiamond,-fZCutDiamond+binsize*nbins); + fGenerForDownScale=new TRandom3(987654321); } //______________________________________________________________________ @@ -84,6 +96,7 @@ AliITSVertexer3D::~AliITSVertexer3D() { // Destructor fLines.Clear("C"); if(fZHisto) delete fZHisto; + if(fGenerForDownScale) delete fGenerForDownScale; } //______________________________________________________________________ @@ -106,10 +119,19 @@ AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree) fCurrentVertex = NULL; Int_t nolines = FindTracklets(itsClusterTree,0); + Int_t rc; if(nolines>=2){ - Int_t rc=Prepare3DVertex(0); - if(fPileupAlgo == 2 && rc == 0) FindVertex3DIterative(); - else if(fPileupAlgo<2 && rc == 0) FindVertex3D(itsClusterTree); + rc=Prepare3DVertex(0); + if(fVert3D.GetNContributors()>0){ + fLines.Clear("C"); + nolines = FindTracklets(itsClusterTree,1); + if(nolines>=2){ + rc=Prepare3DVertex(1); + if(fPileupAlgo == 2 && rc == 0) FindVertex3DIterative(); + else if(fPileupAlgo!=2 && rc == 0) FindVertex3D(itsClusterTree); + if(rc!=0) fVert3D.SetNContributors(0); // exclude this vertex + } + } } if(!fCurrentVertex){ @@ -126,38 +148,21 @@ AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree) Double_t chi2=99999.; Int_t nContr=vtxz->GetNContributors(); fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); + fCurrentVertex->SetDispersion(vtxz->GetDispersion()); fCurrentVertex->SetTitle("vertexer: Z"); fCurrentVertex->SetName("SPDVertexZ"); delete vtxz; } - } - FindMultiplicity(itsClusterTree); + } + if(fComputeMultiplicity) FindMultiplicity(itsClusterTree); return fCurrentVertex; } //______________________________________________________________________ void AliITSVertexer3D::FindVertex3D(TTree *itsClusterTree){ - // 3D algorithm - /* uncomment to debug - printf("Vertex found in first iteration:\n"); - fVert3D.Print(); - printf("Start second iteration\n"); - end of debug lines */ - if(fVert3D.GetNContributors()>0){ - fLines.Clear("C"); - Int_t nolines = FindTracklets(itsClusterTree,1); - if(nolines>=2){ - Int_t rc=Prepare3DVertex(1); - if(rc!=0) fVert3D.SetNContributors(0); // exclude this vertex - } - } - /* uncomment to debug - printf("Vertex found in second iteration:\n"); - fVert3D.Print(); - end of debug lines */ - - Float_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv()); + + Double_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv()); if(vRadius0){ Double_t position[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()}; Double_t covmatrix[6]; @@ -173,6 +178,7 @@ void AliITSVertexer3D::FindVertex3D(TTree *itsClusterTree){ switch(fPileupAlgo){ case 0: PileupFromZ(); break; case 1: FindOther3DVertices(itsClusterTree); break; + case 3: break; // no pileup algo default: AliError("Wrong pileup algorithm"); break; } if(fNoVertices==1){ @@ -184,6 +190,216 @@ void AliITSVertexer3D::FindVertex3D(TTree *itsClusterTree){ //______________________________________________________________________ void AliITSVertexer3D::FindVertex3DIterative(){ + // + + Int_t nLines=fLines.GetEntriesFast(); + Int_t maxPoints=nLines*(nLines-1)/2; + Double_t* xP=new Double_t[maxPoints]; + Double_t* yP=new Double_t[maxPoints]; + Double_t* zP=new Double_t[maxPoints]; + Int_t* index1=new Int_t[maxPoints]; + Int_t* index2=new Int_t[maxPoints]; + Double_t xbeam=fVert3D.GetXv(); + Double_t ybeam=fVert3D.GetYv(); + + Int_t iPoint=0; + for(Int_t ilin1=0; ilin1GetDCA(l2); + if(dca > fDCAcut || dca<0.00001) continue; + Double_t point[3]; + Int_t retc = l1->Cross(l2,point); + if(retc<0)continue; + Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]); + if(rad>fCoarseMaxRCut)continue; + Double_t distFromBeam=TMath::Sqrt((point[0]-xbeam)*(point[0]-xbeam)+(point[1]-ybeam)*(point[1]-ybeam)); + if(distFromBeam>fMaxRCut2) continue; + xP[iPoint]=point[0]; + yP[iPoint]=point[1]; + zP[iPoint]=point[2]; + index1[iPoint]=ilin1; + index2[iPoint]=ilin2; + iPoint++; + } + } + Int_t npoints=iPoint++; + Int_t index=0; + Short_t* mask=new Short_t[npoints]; + for(Int_t ip=0;ip0) cap=isIndUsed[sortedIndex[ind1-1]]; + UInt_t bigger=0; + Int_t biggerindex=-1; + for(Int_t ind2=0;ind2bigger && isIndUsed[ind2]<=cap){ + bigger=isIndUsed[ind2]; + biggerindex=ind2; + } + } + sortedIndex[ind1]=biggerindex; + } + AliDebug(3,Form("Number of clusters before merging = %d\n",nClusters)); + + // Assign lines to clusters/vertices and merge clusters which share 1 line + Int_t nClustersAfterMerge=nClusters; + Int_t* belongsTo=new Int_t[nLines]; + for(Int_t ilin=0; ilin1 associated tracklet) + UInt_t nGoodVert=0; + for(Int_t iclu=0;iclu1) nGoodVert++; + } + + AliDebug(1,Form("Number of good vertices = %d\n",nGoodVert)); + // Calculate vertex coordinates for each cluster + if(nGoodVert>0){ + fVertArray = new AliESDVertex[nGoodVert]; + Int_t iVert=0; + for(Int_t iclu=0;iclu1){ + AliStrLine **arrlin = new AliStrLine*[size]; + Int_t nFilled=0; + for(Int_t ilin=0; ilin(fLines[ilin]); + } + } + AliDebug(3,Form("Vertex %d N associated tracklets = %d out of %d\n",iVert,size,nFilled)); + + fVertArray[iVert]=AliVertexerTracks::TrackletVertexFinder(arrlin,nFilled); + Double_t peak[3]; + fVertArray[iVert].GetXYZ(peak); + AliStrLine **arrlin2 = new AliStrLine*[size]; + Int_t nFilled2=0; + for(Int_t i=0; iGetDistFromPoint(peak)< fDCAcut) + arrlin2[nFilled2++] = dynamic_cast(l1); + } + if(nFilled2>1){ + AliDebug(3,Form("Vertex %d recalculated with %d tracklets\n",iVert,nFilled2)); + fVertArray[iVert]=AliVertexerTracks::TrackletVertexFinder(arrlin2,nFilled2); + } + delete [] arrlin; + delete [] arrlin2; + ++iVert; + } + } + + if(nGoodVert > 1){ + fIsPileup = kTRUE; + fNTrpuv = fVertArray[1].GetNContributors(); + fZpuv = fVertArray[1].GetZv(); + } + + Double_t vRadius=TMath::Sqrt(fVertArray[0].GetXv()*fVertArray[0].GetXv()+fVertArray[0].GetYv()*fVertArray[0].GetYv()); + if(vRadius0){ + Double_t position[3]={fVertArray[0].GetXv(),fVertArray[0].GetYv(),fVertArray[0].GetZv()}; + Double_t covmatrix[6]; + fVertArray[0].GetCovMatrix(covmatrix); + Double_t chi2=99999.; + Int_t nContr=fVertArray[0].GetNContributors(); + fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); + fCurrentVertex->SetTitle("vertexer: 3D"); + fCurrentVertex->SetName("SPDVertex3D"); + fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion()); + } + } + + delete [] index1; + delete [] index2; + delete [] mask; + delete [] isIndUsed; + delete [] sortedIndex; + delete [] belongsTo; + delete [] cluSize; + delete [] xP; + delete [] yP; + delete [] zP; +} +//______________________________________________________________________ +void AliITSVertexer3D::FindVertex3DIterativeMM(){ // Defines the AliESDVertex for the current event Int_t numsor=fLines.GetEntriesFast()*(fLines.GetEntriesFast()-1)/2; //cout<<"AliITSVertexer3D::FindVertexForCurentEvent: Number of tracklets selected for vertexing "<0){ Double_t position[3]={fVertArray[0].GetXv(),fVertArray[0].GetYv(),fVertArray[0].GetZv()}; Double_t covmatrix[6]; @@ -262,101 +478,99 @@ Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){ // considered. Straight lines (=tracklets)are formed. // The tracklets are processed in Prepare3DVertex - if(!GetDetTypeRec())AliFatal("DetTypeRec pointer has not been set"); - - TTree *tR = itsClusterTree; - fDetTypeRec->ResetRecPoints(); - fDetTypeRec->SetTreeAddressR(tR); TClonesArray *itsRec = 0; if(optCuts==0) fZHisto->Reset(); // gc1 are local and global coordinates for layer 1 - Float_t gc1[3]={0.,0.,0.}; + Float_t gc1f[3]={0.,0.,0.}; + Double_t gc1[3]={0.,0.,0.}; // gc2 are local and global coordinates for layer 2 - Float_t gc2[3]={0.,0.,0.}; - - itsRec = fDetTypeRec->RecPoints(); - TBranch *branch = NULL; - branch = tR->GetBranch("ITSRecPoints"); - if(!branch){ - AliError("Null pointer for RecPoints branch"); + Float_t gc2f[3]={0.,0.,0.}; + Double_t gc2[3]={0.,0.,0.}; + AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance(); + itsRec=rpcont->FetchClusters(0,itsClusterTree); + if(!rpcont->IsSPDActive()){ + AliWarning("No SPD rec points found, 3D vertex not calculated"); return -1; } // Set values for cuts - Float_t xbeam=GetNominalPos()[0]; - Float_t ybeam=GetNominalPos()[1]; - Float_t zvert=0.; - Float_t deltaPhi=fCoarseDiffPhiCut; - Float_t deltaR=fCoarseMaxRCut; - Float_t dZmax=fZCutDiamond; - if(fPileupAlgo == 2){ - deltaPhi=fFineDiffPhiCut; - deltaR=fMaxRCut; - if(optCuts != 0)AliWarning(Form("fPileupAlgo=2 AND optCuts=%d has been selected. It should be 0",optCuts)); - } else if(optCuts==1){ + Double_t xbeam=GetNominalPos()[0]; + Double_t ybeam=GetNominalPos()[1]; + Double_t zvert=0.; + Double_t deltaPhi=fCoarseDiffPhiCut; + Double_t deltaR=fCoarseMaxRCut; + Double_t dZmax=fZCutDiamond; + if(optCuts==1){ xbeam=fVert3D.GetXv(); ybeam=fVert3D.GetYv(); zvert=fVert3D.GetZv(); deltaPhi = fDiffPhiMax; deltaR=fMaxRCut; dZmax=fMaxZCut; + if(fPileupAlgo == 2){ + dZmax=fZCutDiamond; + deltaR=fMaxRCut2; + } } else if(optCuts==2){ xbeam=fVert3D.GetXv(); ybeam=fVert3D.GetYv(); - deltaPhi = fDiffPhiMax; + deltaPhi = fDiffPhiforPileup; deltaR=fMaxRCut; } Int_t nrpL1 = 0; // number of rec points on layer 1 Int_t nrpL2 = 0; // number of rec points on layer 2 - - // By default irstL1=0 and lastL1=79 - Int_t firstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1); - Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1; - for(Int_t module= firstL1; module<=lastL1;module++){ // count number of recopints on layer 1 - branch->GetEvent(module); - nrpL1+= itsRec->GetEntries(); - fDetTypeRec->ResetRecPoints(); - } - //By default firstL2=80 and lastL2=239 - Int_t firstL2 = AliITSgeomTGeo::GetModuleIndex(2,1,1); - Int_t lastL2 = AliITSgeomTGeo::GetModuleIndex(3,1,1)-1; - for(Int_t module= firstL2; module<=lastL2;module++){ // count number of recopints on layer 2 - branch->GetEvent(module); - nrpL2+= itsRec->GetEntries(); - fDetTypeRec->ResetRecPoints(); - } + nrpL1=rpcont->GetNClustersInLayerFast(1); + nrpL2=rpcont->GetNClustersInLayerFast(2); if(nrpL1 == 0 || nrpL2 == 0){ + AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",nrpL1,nrpL2)); return -1; } AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2)); + fDoDownScale=kFALSE; + Float_t factDownScal=1.; + Int_t origLaddersOnLayer2=fLadOnLay2; + + if(nrpL1>fMaxNumOfCl || nrpL2>fMaxNumOfCl){ + SetLaddersOnLayer2(2); + fDoDownScale=kTRUE; + factDownScal=(Float_t)fMaxNumOfCl*(Float_t)fMaxNumOfCl/(Float_t)nrpL1/(Float_t)nrpL2; + if(optCuts==1){ + factDownScal*=(fCoarseDiffPhiCut/fDiffPhiMax)*10; + if(factDownScal>1.){ + fDoDownScale=kFALSE; + SetLaddersOnLayer2(origLaddersOnLayer2); + } + } + else if(optCuts==2) return -1; + if(fDoDownScale)AliDebug(1,Form("Too many recpoints on SPD(%d %d ), downscale by %f",nrpL1,nrpL2,factDownScal)); + } Double_t a[3]={xbeam,ybeam,0.}; Double_t b[3]={xbeam,ybeam,10.}; AliStrLine zeta(a,b,kTRUE); - Float_t bField=AliTracker::GetBz()/10.; //T + static Double_t bField=TMath::Abs(AliTracker::GetBz()/10.); //T SetMeanPPtSelTracks(bField); Int_t nolines = 0; // Loop on modules of layer 1 + Int_t firstL1 = TMath::Max(0,AliITSgeomTGeo::GetModuleIndex(1,1,1)); + Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1; for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1 if(!fUseModule[modul1]) continue; - UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1 - branch->GetEvent(modul1); - Int_t nrecp1 = itsRec->GetEntries(); - static TClonesArray prpl1("AliITSRecPoint",nrecp1); - prpl1.SetOwner(); - for(Int_t j=0;jAt(j); - new(prpl1[j])AliITSRecPoint(*recp); - } - fDetTypeRec->ResetRecPoints(); + + UShort_t ladder=modul1/4+1; // ladders are numbered starting from 1 + TClonesArray *prpl1=rpcont->UncheckedGetClusters(modul1); + Int_t nrecp1 = prpl1->GetEntries(); for(Int_t j=0;jkMaxCluPerMod) continue; UShort_t idClu1=modul1*kMaxCluPerMod+j; if(fUsedCluster.TestBitNumber(idClu1)) continue; - AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1.At(j); - recp1->GetGlobalXYZ(gc1); + if(fDoDownScale && fGenerForDownScale->Rndm()>factDownScal) continue; + AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1->At(j); + recp1->GetGlobalXYZ(gc1f); + for(Int_t ico=0;ico<3;ico++)gc1[ico]=gc1f[ico]; + Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam); if(phi1<0)phi1=2*TMath::Pi()+phi1; for(Int_t ladl2=0 ; ladl2AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2); Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1); + if(modul2<0)continue; if(!fUseModule[modul2]) continue; - branch->GetEvent(modul2); + itsRec=rpcont->UncheckedGetClusters(modul2); Int_t nrecp2 = itsRec->GetEntries(); for(Int_t j2=0;j2kMaxCluPerMod) continue; UShort_t idClu2=modul2*kMaxCluPerMod+j2; if(fUsedCluster.TestBitNumber(idClu2)) continue; + AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2); - recp2->GetGlobalXYZ(gc2); + recp2->GetGlobalXYZ(gc2f); + for(Int_t ico=0;ico<3;ico++)gc2[ico]=gc2f[ico]; Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam); if(phi2<0)phi2=2*TMath::Pi()+phi2; Double_t diff = TMath::Abs(phi2-phi1); if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; - if(optCuts==0 && diffFill(zr0); } if(diff>deltaPhi)continue; @@ -396,6 +613,7 @@ Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){ Double_t deltaZ=cp[2]-zvert; if(TMath::Abs(deltaZ)>dZmax)continue; + if(nolines == 0){ if(fLines.GetEntriesFast()>0)fLines.Clear("C"); } @@ -403,50 +621,54 @@ Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){ recp2->GetGlobalCov(cov); - Float_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); - Float_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); - Float_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction + Double_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); + Double_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); + Double_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction - Float_t curvErr=0; + Double_t curvErr=0; if(bField>0.00001){ - Float_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm - Float_t dRad=TMath::Sqrt(TMath::Power((gc1[0]-gc2[0]),2)+TMath::Power((gc1[1]-gc2[1]),2)); - Float_t aux=dRad/2.+rad1; + Double_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm + Double_t dRad=TMath::Sqrt((gc1[0]-gc2[0])*(gc1[0]-gc2[0])+(gc1[1]-gc2[1])*(gc1[1]-gc2[1])); + Double_t aux=dRad/2.+rad1; curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm } - Float_t sigmasq[3]; + Double_t sigmasq[3]; sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor; sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor; sigmasq[2]=cov[5]*factor*factor; // Multiple scattering - Float_t pOverMass=fMeanPSelTrk/0.140; - Float_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass); - Float_t p2=fMeanPSelTrk*fMeanPSelTrk; - Float_t rBP=GetPipeRadius(); - Float_t dBP=0.08/35.3; // 800 um of Be - Float_t dL1=0.01; //approx. 1% of radiation length - Float_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP; - Float_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1; - Float_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1)); - Float_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP)); + Double_t pOverMass=fMeanPSelTrk/0.140; + Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass); + Double_t p2=fMeanPSelTrk*fMeanPSelTrk; + Double_t rBP=GetPipeRadius(); + 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++){ sigmasq[ico]+=rtantheta1*rtantheta1*factor*factor/3.; sigmasq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.; } - Float_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.}; + Double_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.}; if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0]; if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1]; if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2]; new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE,idClu1,idClu2); } - fDetTypeRec->ResetRecPoints(); } } } - prpl1.Clear(); } + + if(fDoDownScale){ + SetLaddersOnLayer2(origLaddersOnLayer2); + } + + if(nolines == 0)return -2; return nolines; } @@ -456,42 +678,46 @@ Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ // Finds the 3D vertex information using tracklets Int_t retcode = -1; - Float_t xbeam=GetNominalPos()[0]; - Float_t ybeam=GetNominalPos()[1]; - Float_t zvert=0.; - Float_t deltaR=fCoarseMaxRCut; - if(fPileupAlgo == 2) { - deltaR=fMaxRCut; - if(optCuts!=0)AliWarning(Form("fPileupAlgo=2 AND optCuts=%d. It should be 0",optCuts)); - } - Float_t dZmax=fZCutDiamond; + Double_t xbeam=GetNominalPos()[0]; + Double_t ybeam=GetNominalPos()[1]; + Double_t zvert=0.; + Double_t deltaR=fCoarseMaxRCut; + Double_t dZmax=fZCutDiamond; if(optCuts==1){ xbeam=fVert3D.GetXv(); ybeam=fVert3D.GetYv(); zvert=fVert3D.GetZv(); deltaR=fMaxRCut; dZmax=fMaxZCut; + if(fPileupAlgo == 2){ + dZmax=fZCutDiamond; + deltaR=fMaxRCut2; + } }else if(optCuts==2){ xbeam=fVert3D.GetXv(); ybeam=fVert3D.GetYv(); deltaR=fMaxRCut; } - Float_t rl=-fCoarseMaxRCut; - Float_t rh=fCoarseMaxRCut; - Float_t zl=-fZCutDiamond; - Float_t zh=fZCutDiamond; + Double_t origBinSizeR=fBinSizeR; + Double_t origBinSizeZ=fBinSizeZ; + if(fDoDownScale){ + SetBinSizeR(0.05); + SetBinSizeZ(0.05); + } + + Double_t rl=-fCoarseMaxRCut; + Double_t rh=fCoarseMaxRCut; + Double_t zl=-fZCutDiamond; + Double_t zh=fZCutDiamond; Int_t nbr=(Int_t)((rh-rl)/fBinSizeR+0.0001); Int_t nbz=(Int_t)((zh-zl)/fBinSizeZ+0.0001); Int_t nbrcs=(Int_t)((rh-rl)/(fBinSizeR*2.)+0.0001); Int_t nbzcs=(Int_t)((zh-zl)/(fBinSizeZ*2.)+0.0001); - TH3F *h3d = NULL; - TH3F *h3dcs = NULL; - if(fPileupAlgo !=2){ - h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh); - h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh); - } + TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh); + TH3F *h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh); + // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed) Int_t *validate = new Int_t [fLines.GetEntriesFast()]; for(Int_t i=0; iGetDCA(l2); + Double_t dca=l1->GetDCA(l2); if(dca > fDCAcut || dca<0.00001) continue; Double_t point[3]; Int_t retc = l1->Cross(l2,point); @@ -514,10 +740,8 @@ Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ if(raddist>deltaR)continue; validate[i]=1; validate[j]=1; - if(fPileupAlgo != 2){ - h3d->Fill(point[0],point[1],point[2]); - h3dcs->Fill(point[0],point[1],point[2]); - } + h3d->Fill(point[0],point[1],point[2]); + h3dcs->Fill(point[0],point[1],point[2]); } } @@ -527,10 +751,8 @@ Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ for(Int_t i=0; i=1)numbtracklets++; if(numbtracklets<2){ delete [] validate; - if(fPileupAlgo != 2){ - delete h3d; - delete h3dcs; - } + delete h3d; + delete h3dcs; return retcode; } @@ -541,10 +763,12 @@ Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast())); delete [] validate; - // Exit here if Pileup Algorithm 2 has been chosen - if(fPileupAlgo == 2)return 0; - // - + // Exit here if Pileup Algorithm 2 has been chosen during second loop + if(fPileupAlgo == 2 && optCuts==1){ + delete h3d; + delete h3dcs; + return 0; + } // Find peaks in histos @@ -552,8 +776,8 @@ Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ Int_t ntrkl,ntimes; FindPeaks(h3d,peak,ntrkl,ntimes); delete h3d; - Float_t binsizer=(rh-rl)/nbr; - Float_t binsizez=(zh-zl)/nbz; + Double_t binsizer=(rh-rl)/nbr; + Double_t binsizez=(zh-zl)/nbz; if(optCuts==0 && (ntrkl<=2 || ntimes>1)){ ntrkl=0; ntimes=0; @@ -564,9 +788,50 @@ Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ } delete h3dcs; + // Finer Histo in limited range in case of high mult. + if(fDoDownScale){ + SetBinSizeR(0.01); + SetBinSizeZ(0.01); + Double_t xl=peak[0]-0.3; + Double_t xh=peak[0]+0.3; + Double_t yl=peak[1]-0.3; + Double_t yh=peak[1]+0.3; + zl=peak[2]-0.5; + zh=peak[2]+0.5; + Int_t nbxfs=(Int_t)((xh-xl)/fBinSizeR+0.0001); + Int_t nbyfs=(Int_t)((yh-yl)/fBinSizeR+0.0001); + Int_t nbzfs=(Int_t)((zh-zl)/fBinSizeZ+0.0001); + + TH3F *h3dfs = new TH3F("h3dfs","xyz distribution",nbxfs,xl,xh,nbyfs,yl,yh,nbzfs,zl,zh); + for(Int_t i=0; iCross(l2,point); + if(retc<0)continue; + h3dfs->Fill(point[0],point[1],point[2]); + } + } + ntrkl=0; + ntimes=0; + + Double_t newpeak[3]={0.,0.,0.}; + FindPeaks(h3dfs,newpeak,ntrkl,ntimes); + if(ntimes==1){ + for(Int_t iCoo=0; iCoo<3; iCoo++) peak[iCoo]=newpeak[iCoo]; + binsizer=fBinSizeR; + binsizez=fBinSizeZ; + } + delete h3dfs; + SetBinSizeR(origBinSizeR); + SetBinSizeZ(origBinSizeZ); + } + + // Second selection loop - Float_t bs=(binsizer+binsizez)/2.; + Double_t bs=(binsizer+binsizez)/2.; for(Int_t i=0; iGetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i); @@ -581,10 +846,33 @@ Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ // make a further selection on tracklets based on this first candidate fVert3D.GetXYZ(peak); AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2])); + Int_t *validate2 = new Int_t [fLines.GetEntriesFast()]; + for(Int_t i=0; iGetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i); + if(optCuts==2){ // temporarily only for pileup + for(Int_t j=i+1; jGetDCA(l2)<0.00001){ + Int_t delta1=(Int_t)l1->GetIdPoint(0)-(Int_t)l2->GetIdPoint(0); + Int_t delta2=(Int_t)l1->GetIdPoint(1)-(Int_t)l2->GetIdPoint(1); + Int_t deltamod1=(Int_t)l1->GetIdPoint(0)/kMaxCluPerMod + -(Int_t)l2->GetIdPoint(0)/kMaxCluPerMod; + Int_t deltamod2=(Int_t)l1->GetIdPoint(1)/kMaxCluPerMod + -(Int_t)l2->GetIdPoint(1)/kMaxCluPerMod; + // remove tracklets sharing a point + if( (delta1==0 && deltamod2==0) || + (delta2==0 && deltamod1==0) ) validate2[j]=0; + } + } + } } + for(Int_t i=0; i1){// this new tracklet selection is used @@ -595,7 +883,7 @@ Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ } //________________________________________________________ -void AliITSVertexer3D::SetMeanPPtSelTracks(Float_t fieldTesla){ +void AliITSVertexer3D::SetMeanPPtSelTracks(Double_t fieldTesla){ // Sets mean values of Pt based on the field // for P (used in multiple scattering) the most probable value is used if(TMath::Abs(fieldTesla-0.5)<0.01){ @@ -631,11 +919,11 @@ void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklet Int_t peak2bin[3]={-1,-1,-1}; Int_t bc2=-1; for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){ - Float_t xval = xax->GetBinCenter(i); + Double_t xval = xax->GetBinCenter(i); for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){ - Float_t yval = yax->GetBinCenter(j); + Double_t yval = yax->GetBinCenter(j); for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){ - Float_t zval = zax->GetBinCenter(k); + Double_t zval = zax->GetBinCenter(k); Int_t bc =(Int_t)histo->GetBinContent(i,j,k); if(bc==0) continue; if(bc>nOfTracklets){ @@ -736,7 +1024,7 @@ void AliITSVertexer3D::PileupFromZ(){ Int_t nPeaks=AliITSVertexerZ::GetPeakRegion(fZHisto,binmin,binmax); if(nPeaks==2)AliWarning("2 peaks found"); Int_t firstPeakCont=0; - Float_t firstPeakPos=0.; + Double_t firstPeakPos=0.; for(Int_t i=binmin-1;i<=binmax+1;i++){ firstPeakCont+=static_cast(fZHisto->GetBinContent(i)); firstPeakPos+=fZHisto->GetBinContent(i)*fZHisto->GetBinCenter(i); @@ -763,14 +1051,21 @@ void AliITSVertexer3D::PileupFromZ(){ //________________________________________________________ void AliITSVertexer3D::PrintStatus() const { // Print current status - printf("=======================================================\n"); - printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut); - printf("Cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut); - printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut); + printf("========= First step selections =====================\n"); printf("Cut on diamond (Z) %f\n",fZCutDiamond); + printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut); + printf("Loose cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut); printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut); + printf("========= Second step selections ====================\n"); + printf("Cut on tracklet-to-first-vertex Z distance %f\n",fMaxZCut); printf("Max Phi difference: %f\n",fDiffPhiMax); + printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut); + printf("Cut on tracklet DCA to beam axis (algo2) %f\n",fMaxRCut2); + printf("========= Pileup selections =========================\n"); printf("Pileup algo: %d\n",fPileupAlgo); - printf("Min DCA to 1st vetrtex for pileup: %f\n",fDCAforPileup); + printf("Min DCA to 1st vertex for pileup (algo 0 and 1): %f\n",fDCAforPileup); + printf("Cut on distance between pair-vertices (algo 2): %f\n",fCutOnPairs); + printf("Maximum number of clusters allowed on L1 or L2: %d\n",fMaxNumOfCl); printf("=======================================================\n"); } +