X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSVertexerZ.cxx;h=798b9fb1b3ed5cc65ddd1a4df9ccc430085d26d5;hb=f5e33cc07e8750f569fba3903ce6f387d30d46bf;hp=dcb695d6458892ba54ee130a4dd1dc572f08ed04;hpb=8c32ba4472d30b86eeaf84510539ce57bdd5d176;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSVertexerZ.cxx b/ITS/AliITSVertexerZ.cxx index dcb695d6458..798b9fb1b3e 100644 --- a/ITS/AliITSVertexerZ.cxx +++ b/ITS/AliITSVertexerZ.cxx @@ -20,8 +20,8 @@ #include #include "AliESDVertex.h" #include "AliITSgeomTGeo.h" -#include "AliITSDetTypeRec.h" #include "AliITSRecPoint.h" +#include "AliITSRecPointContainer.h" #include "AliITSZPoint.h" ///////////////////////////////////////////////////////////////// @@ -32,6 +32,8 @@ // It can be used successfully with Pb-Pb collisions //////////////////////////////////////////////////////////////// +using std::endl; +using std::cout; ClassImp(AliITSVertexerZ) @@ -51,7 +53,9 @@ fHighLim(0.), fStepCoarse(0), fTolerance(0.), fMaxIter(0), -fWindowWidth(0) { +fWindowWidth(0), +fSearchForPileup(kTRUE) +{ // Default constructor SetDiffPhiMax(); SetFirstLayerModules(); @@ -80,7 +84,9 @@ fHighLim(0.), fStepCoarse(0), fTolerance(0.), fMaxIter(0), -fWindowWidth(0) { +fWindowWidth(0), +fSearchForPileup(kTRUE) +{ // Standard Constructor SetDiffPhiMax(); SetFirstLayerModules(); @@ -106,7 +112,7 @@ AliITSVertexerZ::~AliITSVertexerZ() { void AliITSVertexerZ::ConfigIterations(Int_t noiter,Float_t *ptr){ // configure the iterative procedure to gain efficiency for // pp events with very low multiplicity - Float_t defaults[5]={0.05,0.1,0.2,0.3,0.5}; + Float_t defaults[5]={0.02,0.05,0.1,0.2,0.3}; fMaxIter=noiter; if(noiter>5){ Error("ConfigIterations","Maximum number of iterations is 5\n"); @@ -180,47 +186,42 @@ AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(TTree *itsClusterTree){ fDiffPhiMax=diffPhiMaxOrig; } } - FindMultiplicity(itsClusterTree); + if(fComputeMultiplicity) FindMultiplicity(itsClusterTree); return fCurrentVertex; } //______________________________________________________________________ void AliITSVertexerZ::VertexZFinder(TTree *itsClusterTree){ // Defines the AliESDVertex for the current event + delete fCurrentVertex; fCurrentVertex = 0; - - TTree *tR = itsClusterTree; - fDetTypeRec->SetTreeAddressR(tR); + Double_t startPos[3]={GetNominalPos()[0],GetNominalPos()[1],GetNominalPos()[2]}; + Double_t startCov[6]={GetNominalCov()[0],GetNominalCov()[1],GetNominalCov()[2], + GetNominalCov()[3],GetNominalCov()[4],GetNominalCov()[5]}; + ResetVertex(); TClonesArray *itsRec = 0; // lc1 and gc1 are local and global coordinates for layer 1 - Float_t lc1[3]; for(Int_t ii=0; ii<3; ii++) lc1[ii]=0.; - Float_t gc1[3]; for(Int_t ii=0; ii<3; ii++) gc1[ii]=0.; + Float_t gc1[3]={0.,0.,0.}; // ; for(Int_t ii=0; ii<3; ii++) gc1[ii]=0.; // lc2 and gc2 are local and global coordinates for layer 2 - Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.; - Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.; - - itsRec = fDetTypeRec->RecPoints(); - TBranch *branch; - branch = tR->GetBranch("ITSRecPoints"); + Float_t gc2[3]={0.,0.,0.}; //; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.; + AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance(); + rpcont->FetchClusters(0,itsClusterTree); + if(!rpcont->IsSPDActive()){ + AliWarning("Null pointer for RecPoints branch, vertex not calculated"); + ResetHistograms(); + fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-2); + return; + } Int_t nrpL1 = 0; Int_t nrpL2 = 0; + nrpL1=rpcont->GetNClustersInLayerFast(1); + nrpL2=rpcont->GetNClustersInLayerFast(2); - // By default fFirstL1=0 and fLastL1=79 - for(Int_t module= fFirstL1; module<=fLastL1;module++){ - branch->GetEvent(module); - nrpL1+= itsRec->GetEntries(); - fDetTypeRec->ResetRecPoints(); - } - //By default fFirstL2=80 and fLastL2=239 - for(Int_t module= fFirstL2; module<=fLastL2;module++){ - branch->GetEvent(module); - nrpL2+= itsRec->GetEntries(); - fDetTypeRec->ResetRecPoints(); - } if(nrpL1 == 0 || nrpL2 == 0){ + AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",nrpL1,nrpL2)); ResetHistograms(); - fCurrentVertex = new AliESDVertex(0.,5.3,-2); + fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-2); return; } // Force a coarse bin size of 200 microns if the number of clusters on layer 2 @@ -242,81 +243,74 @@ void AliITSVertexerZ::VertexZFinder(TTree *itsClusterTree){ Float_t theta2BP=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dBP); Float_t theta2L1=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dL1); */ + Int_t nEntriesMod[kNSPDMod]; + TClonesArray* recpArr[kNSPDMod]; + for(Int_t modul=0; modulUncheckedGetClusters(modul); + nEntriesMod[modul]=recpArr[modul]->GetEntriesFast(); + } + } + Int_t maxdim=TMath::Min(nrpL1*nrpL2,50000); // temporary; to limit the size in PbPb static TClonesArray points("AliITSZPoint",maxdim); Int_t nopoints =0; for(Int_t modul1= fFirstL1; modul1<=fLastL1;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(); + TClonesArray *prpl1=recpArr[modul1]; //rpcont->UncheckedGetClusters(modul1); + Int_t nrecp1 = nEntriesMod[modul1]; //prpl1->GetEntries(); for(Int_t j1=0;j1GetDetLocalX(); - lc1[2]=recp->GetDetLocalZ(); - geom->LtoG(modul1,lc1,gc1); - // Global coordinates of this recpoints - */ - recp->GetGlobalXYZ(gc1); + AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1->At(j1); + recp1->GetGlobalXYZ(gc1); gc1[0]-=GetNominalPos()[0]; // Possible beam offset in the bending plane gc1[1]-=GetNominalPos()[1]; // " " - Float_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); Float_t phi1 = TMath::ATan2(gc1[1],gc1[0]); - if(phi1<0)phi1+=2*TMath::Pi(); - Float_t zc1=gc1[2]; - Float_t erz1=recp->GetSigmaZ2(); + if(phi1<0)phi1+=TMath::TwoPi(); 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); - Int_t nrecp2 = itsRec->GetEntries(); + itsRec=recpArr[modul2]; // rpcont->UncheckedGetClusters(modul2); + Int_t nrecp2 = nEntriesMod[modul2]; // itsRec->GetEntries(); for(Int_t j2=0;j2At(j2); - /* - lc2[0]=recp->GetDetLocalX(); - lc2[2]=recp->GetDetLocalZ(); - geom->LtoG(modul2,lc2,gc2); - */ - recp->GetGlobalXYZ(gc2); + AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2); + recp2->GetGlobalXYZ(gc2); gc2[0]-=GetNominalPos()[0]; gc2[1]-=GetNominalPos()[1]; - Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]); - if(phi2<0)phi2+=2*TMath::Pi(); - Float_t zc2=gc2[2]; - Float_t erz2=recp->GetSigmaZ2(); + if(phi2<0)phi2+=TMath::TwoPi(); Float_t diff = TMath::Abs(phi2-phi1); - if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; + if(diff>TMath::Pi())diff=TMath::TwoPi()-diff; if(diffGetSigmaZ2(); + Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); + Float_t zc2=gc2[2]; + Float_t erz2=recp2->GetSigmaZ2(); // Float_t tgth=(zc2[j]-zc1[i])/(r2-r1); // slope (used for multiple scattering) Float_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius - Float_t ezr0q=(r2*r2*erz1+r1*r1*erz2)/(r2-r1)/(r2-r1); //error on Z @ null radius - /* - // Multiple scattering - ezr0q+=r1*r1*(1+tgth*tgth)*theta2L1/2; // multiple scattering in layer 1 - ezr0q+=rBP*rBP*(1+tgth*tgth)*theta2BP/2; // multiple scattering in beam pipe - */ + Float_t ezr0q=(r2*r2*erz1+r1*r1*erz2)/((r2-r1)*(r2-r1)); //error on Z @ null radius + /* + // Multiple scattering + ezr0q+=r1*r1*(1+tgth*tgth)*theta2L1/2; // multiple scattering in layer 1 + ezr0q+=rBP*rBP*(1+tgth*tgth)*theta2BP/2; // multiple scattering in beam pipe + */ if(nopointsFill(zr0); } } - fDetTypeRec->ResetRecPoints(); } } } - prpl1.Clear(); } points.Sort(); @@ -325,7 +319,7 @@ void AliITSVertexerZ::VertexZFinder(TTree *itsClusterTree){ if(contents<1.){ // Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n"); ResetHistograms(); - fCurrentVertex = new AliESDVertex(0.,5.3,-1); + fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-1); points.Clear(); return; } @@ -336,11 +330,12 @@ void AliITSVertexerZ::VertexZFinder(TTree *itsClusterTree){ if(hc->GetBinContent(hc->GetMaximumBin())<3)hc->Rebin(4); Int_t binmin,binmax; Int_t nPeaks=GetPeakRegion(hc,binmin,binmax); - if(nPeaks==2)AliWarning("2 peaks found"); + if(nPeaks==2)AliDebug(2,"2 peaks found"); Float_t zm =0.; Float_t ezm =0.; Float_t lim1 = hc->GetBinLowEdge(binmin); Float_t lim2 = hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax); + Float_t widthSR=lim2-lim1; if(nPeaks ==1 && (lim2-lim1)GetZ()>lim1 && p->GetZ()GetErrZ(); @@ -374,36 +369,51 @@ void AliITSVertexerZ::VertexZFinder(TTree *itsClusterTree){ } niter++; } while(niter<10 && TMath::Abs((zm-lim1)-(lim2-zm))>fTolerance); - fCurrentVertex = new AliESDVertex(zm,ezm,ncontr); - fCurrentVertex->SetTitle("vertexer: B"); + if(nPeaks==2) ezm=widthSR; + Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],zm}; + Double_t covmatrix[6]={GetNominalCov()[0],0.,GetNominalCov()[2],0.,0.,ezm}; + fCurrentVertex = new AliESDVertex(position,covmatrix,99999.,ncontr); + fCurrentVertex->SetTitle("vertexer: Z"); + fCurrentVertex->SetDispersion(fDiffPhiMax); + fNoVertices=1; points.Clear(); - fIsPileup=kFALSE; - fNTrpuv=-2; - if(ncontr>fMinTrackletsForPilup){ + if(fSearchForPileup && ncontr>fMinTrackletsForPilup){ Float_t secPeakPos; Int_t ncontr2=FindSecondPeak(fZCombc,binmin,binmax,secPeakPos); if(ncontr2>=fMinTrackletsForPilup){ fIsPileup=kTRUE; + fNoVertices=2; fZpuv=secPeakPos; fNTrpuv=ncontr2; + AliESDVertex secondVert(secPeakPos,0.1,ncontr2); + fVertArray = new AliESDVertex[2]; + fVertArray[0]=(*fCurrentVertex); + fVertArray[1]=secondVert; } - } + } + if(fNoVertices==1){ + fVertArray = new AliESDVertex[1]; + fVertArray[0]=(*fCurrentVertex); + } + ResetHistograms(); return; } //_____________________________________________________________________ -Int_t AliITSVertexerZ::FindSecondPeak(TH1F* h, Int_t binmin,Int_t binmax, Float_t& secPeakPos){ +Int_t AliITSVertexerZ::FindSecondPeak(TH1F* h, Int_t binmin,Int_t binmax, Float_t& secPeakPos){ + // Resets bin contents between binmin and binmax and then search + // for a second peak position for(Int_t i=binmin-1;i<=binmax+1;i++){ h->SetBinContent(i,0.); } Int_t secPeakBin=h->GetMaximumBin(); secPeakPos=h->GetBinCenter(secPeakBin); - Int_t secPeakCont=h->GetBinContent(secPeakBin); - secPeakCont+=h->GetBinContent(secPeakBin-1); - secPeakCont+=h->GetBinContent(secPeakBin+1); - secPeakCont+=h->GetBinContent(secPeakBin-2); - secPeakCont+=h->GetBinContent(secPeakBin+2); + Int_t secPeakCont=(Int_t)h->GetBinContent(secPeakBin); + secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-1); + secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+1); + secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-2); + secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+2); return secPeakCont; }