]>
Commit | Line | Data |
---|---|---|
70c95f95 | 1 | /************************************************************************** |
2 | * Copyright(c) 2006-2008, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
70c95f95 | 15 | #include <TTree.h> |
b21c1af0 | 16 | #include "AliRunLoader.h" |
6a0d56b8 | 17 | #include "AliESDVertex.h" |
cd706e57 | 18 | #include "AliLog.h" |
6a0d56b8 | 19 | #include "AliStrLine.h" |
20 | #include "AliTracker.h" | |
70c95f95 | 21 | #include "AliITSDetTypeRec.h" |
22 | #include "AliITSRecPoint.h" | |
b21c1af0 | 23 | #include "AliITSRecPointContainer.h" |
32e63e47 | 24 | #include "AliITSgeomTGeo.h" |
6a0d56b8 | 25 | #include "AliVertexerTracks.h" |
26 | #include "AliITSVertexer3D.h" | |
0b7bac55 | 27 | #include "AliITSVertexerZ.h" |
6b4d9537 | 28 | #include "AliITSSortTrkl.h" |
70c95f95 | 29 | ///////////////////////////////////////////////////////////////// |
30 | // this class implements a method to determine | |
31 | // the 3 coordinates of the primary vertex | |
9801e044 | 32 | // optimized for |
33 | // p-p collisions | |
70c95f95 | 34 | //////////////////////////////////////////////////////////////// |
35 | ||
9801e044 | 36 | const Int_t AliITSVertexer3D::fgkMaxNumOfClDefault = 500; |
37 | ||
70c95f95 | 38 | ClassImp(AliITSVertexer3D) |
39 | ||
6a0d56b8 | 40 | /* $Id$ */ |
41 | ||
70c95f95 | 42 | //______________________________________________________________________ |
7d766abb | 43 | AliITSVertexer3D::AliITSVertexer3D(): |
44 | AliITSVertexer(), | |
45 | fLines("AliStrLine",1000), | |
46 | fVert3D(), | |
47 | fCoarseDiffPhiCut(0.), | |
48 | fFineDiffPhiCut(0.), | |
49 | fCutOnPairs(0.), | |
50 | fCoarseMaxRCut(0.), | |
51 | fMaxRCut(0.), | |
87104b06 | 52 | fMaxRCut2(0.), |
7d766abb | 53 | fZCutDiamond(0.), |
54 | fMaxZCut(0.), | |
55 | fDCAcut(0.), | |
56 | fDiffPhiMax(0.), | |
57 | fMeanPSelTrk(0.), | |
58 | fMeanPtSelTrk(0.), | |
59 | fUsedCluster(kMaxCluPerMod*kNSPDMod), | |
60 | fZHisto(0), | |
61 | fDCAforPileup(0.), | |
c238ba33 | 62 | fDiffPhiforPileup(0.), |
7d766abb | 63 | fBinSizeR(0.), |
64 | fBinSizeZ(0.), | |
9801e044 | 65 | fPileupAlgo(0), |
66 | fMaxNumOfCl(fgkMaxNumOfClDefault) | |
6a0d56b8 | 67 | { |
70c95f95 | 68 | // Default constructor |
69 | SetCoarseDiffPhiCut(); | |
6b4d9537 | 70 | SetFineDiffPhiCut(); |
71 | SetCutOnPairs(); | |
05d1294c | 72 | SetCoarseMaxRCut(); |
70c95f95 | 73 | SetMaxRCut(); |
87104b06 | 74 | SetMaxRCutAlgo2(); |
70c95f95 | 75 | SetZCutDiamond(); |
05d1294c | 76 | SetMaxZCut(); |
7203e11a | 77 | SetDCACut(); |
70c95f95 | 78 | SetDiffPhiMax(); |
6a0d56b8 | 79 | SetMeanPSelTracks(); |
80 | SetMeanPtSelTracks(); | |
cd1c2af1 | 81 | SetMinDCAforPileup(); |
c238ba33 | 82 | SetDeltaPhiforPileup(); |
cd1c2af1 | 83 | SetPileupAlgo(); |
7d766abb | 84 | SetBinSizeR(); |
85 | SetBinSizeZ(); | |
242c8025 | 86 | Double_t binsize=0.02; // default 200 micron |
cd1c2af1 | 87 | Int_t nbins=static_cast<Int_t>(1+2*fZCutDiamond/binsize); |
88 | fZHisto=new TH1F("hz","",nbins,-fZCutDiamond,-fZCutDiamond+binsize*nbins); | |
70c95f95 | 89 | } |
90 | ||
70c95f95 | 91 | //______________________________________________________________________ |
92 | AliITSVertexer3D::~AliITSVertexer3D() { | |
93 | // Destructor | |
f5f6da22 | 94 | fLines.Clear("C"); |
cd1c2af1 | 95 | if(fZHisto) delete fZHisto; |
70c95f95 | 96 | } |
97 | ||
27167524 | 98 | //______________________________________________________________________ |
99 | void AliITSVertexer3D::ResetVert3D(){ | |
100 | // | |
6b4d9537 | 101 | ResetVertex(); |
27167524 | 102 | fVert3D.SetXv(0.); |
103 | fVert3D.SetYv(0.); | |
104 | fVert3D.SetZv(0.); | |
105 | fVert3D.SetDispersion(0.); | |
106 | fVert3D.SetNContributors(0); | |
cd1c2af1 | 107 | fUsedCluster.ResetAllBits(0); |
27167524 | 108 | } |
70c95f95 | 109 | //______________________________________________________________________ |
308c2f7c | 110 | AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree){ |
70c95f95 | 111 | // Defines the AliESDVertex for the current event |
27167524 | 112 | ResetVert3D(); |
308c2f7c | 113 | AliDebug(1,"FindVertexForCurrentEvent - 3D - PROCESSING NEXT EVENT"); |
6b4d9537 | 114 | fLines.Clear("C"); |
115 | fCurrentVertex = NULL; | |
70c95f95 | 116 | |
308c2f7c | 117 | Int_t nolines = FindTracklets(itsClusterTree,0); |
87104b06 | 118 | Int_t rc; |
0b7bac55 | 119 | if(nolines>=2){ |
87104b06 | 120 | rc=Prepare3DVertex(0); |
121 | if(fVert3D.GetNContributors()>0){ | |
122 | fLines.Clear("C"); | |
46004624 | 123 | nolines = FindTracklets(itsClusterTree,1); |
87104b06 | 124 | if(nolines>=2){ |
125 | rc=Prepare3DVertex(1); | |
126 | if(fPileupAlgo == 2 && rc == 0) FindVertex3DIterative(); | |
127 | else if(fPileupAlgo<2 && rc == 0) FindVertex3D(itsClusterTree); | |
128 | if(rc!=0) fVert3D.SetNContributors(0); // exclude this vertex | |
129 | } | |
130 | } | |
0b7bac55 | 131 | } |
6b4d9537 | 132 | |
0b7bac55 | 133 | if(!fCurrentVertex){ |
134 | AliITSVertexerZ vertz(GetNominalPos()[0],GetNominalPos()[1]); | |
b96ee725 | 135 | vertz.SetDetTypeRec(GetDetTypeRec()); |
0b7bac55 | 136 | AliDebug(1,"Call Vertexer Z\n"); |
cc282e54 | 137 | vertz.SetLowLimit(-fZCutDiamond); |
138 | vertz.SetHighLimit(fZCutDiamond); | |
0b7bac55 | 139 | AliESDVertex* vtxz = vertz.FindVertexForCurrentEvent(itsClusterTree); |
140 | if(vtxz){ | |
141 | Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],vtxz->GetZv()}; | |
142 | Double_t covmatrix[6]; | |
143 | vtxz->GetCovMatrix(covmatrix); | |
144 | Double_t chi2=99999.; | |
145 | Int_t nContr=vtxz->GetNContributors(); | |
146 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
bf2e0ad4 | 147 | fCurrentVertex->SetDispersion(vtxz->GetDispersion()); |
0b7bac55 | 148 | fCurrentVertex->SetTitle("vertexer: Z"); |
149 | fCurrentVertex->SetName("SPDVertexZ"); | |
150 | delete vtxz; | |
0b7bac55 | 151 | } |
152 | ||
70c95f95 | 153 | } |
308c2f7c | 154 | FindMultiplicity(itsClusterTree); |
70c95f95 | 155 | return fCurrentVertex; |
156 | } | |
157 | ||
6b4d9537 | 158 | //______________________________________________________________________ |
159 | void AliITSVertexer3D::FindVertex3D(TTree *itsClusterTree){ | |
87104b06 | 160 | |
0496a544 | 161 | Double_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv()); |
6b4d9537 | 162 | if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){ |
163 | Double_t position[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()}; | |
164 | Double_t covmatrix[6]; | |
165 | fVert3D.GetCovMatrix(covmatrix); | |
166 | Double_t chi2=99999.; | |
167 | Int_t nContr=fVert3D.GetNContributors(); | |
168 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
169 | fCurrentVertex->SetTitle("vertexer: 3D"); | |
170 | fCurrentVertex->SetName("SPDVertex3D"); | |
171 | fCurrentVertex->SetDispersion(fVert3D.GetDispersion()); | |
172 | fNoVertices=1; | |
173 | ||
174 | switch(fPileupAlgo){ | |
175 | case 0: PileupFromZ(); break; | |
176 | case 1: FindOther3DVertices(itsClusterTree); break; | |
177 | default: AliError("Wrong pileup algorithm"); break; | |
178 | } | |
179 | if(fNoVertices==1){ | |
180 | fVertArray = new AliESDVertex[1]; | |
181 | fVertArray[0]=(*fCurrentVertex); | |
182 | } | |
183 | } | |
184 | } | |
185 | ||
186 | //______________________________________________________________________ | |
187 | void AliITSVertexer3D::FindVertex3DIterative(){ | |
87104b06 | 188 | // |
189 | ||
190 | Int_t nLines=fLines.GetEntriesFast(); | |
191 | Int_t maxPoints=nLines*(nLines-1)/2; | |
192 | Double_t* xP=new Double_t[maxPoints]; | |
193 | Double_t* yP=new Double_t[maxPoints]; | |
194 | Double_t* zP=new Double_t[maxPoints]; | |
195 | Int_t* index1=new Int_t[maxPoints]; | |
196 | Int_t* index2=new Int_t[maxPoints]; | |
197 | Double_t xbeam=fVert3D.GetXv(); | |
198 | Double_t ybeam=fVert3D.GetYv(); | |
199 | ||
200 | Int_t iPoint=0; | |
201 | for(Int_t ilin1=0; ilin1<nLines; ilin1++){ | |
202 | AliStrLine *l1 = (AliStrLine*)fLines.At(ilin1); | |
203 | for(Int_t ilin2=ilin1+1; ilin2<nLines; ilin2++){ | |
204 | AliStrLine *l2 = (AliStrLine*)fLines.At(ilin2); | |
205 | Double_t dca=l1->GetDCA(l2); | |
206 | if(dca > fDCAcut || dca<0.00001) continue; | |
207 | Double_t point[3]; | |
208 | Int_t retc = l1->Cross(l2,point); | |
209 | if(retc<0)continue; | |
210 | Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]); | |
211 | if(rad>fCoarseMaxRCut)continue; | |
212 | Double_t distFromBeam=TMath::Sqrt((point[0]-xbeam)*(point[0]-xbeam)+(point[1]-ybeam)*(point[1]-ybeam)); | |
213 | if(distFromBeam>fMaxRCut2) continue; | |
214 | xP[iPoint]=point[0]; | |
215 | yP[iPoint]=point[1]; | |
216 | zP[iPoint]=point[2]; | |
217 | index1[iPoint]=ilin1; | |
218 | index2[iPoint]=ilin2; | |
219 | iPoint++; | |
220 | } | |
221 | } | |
222 | Int_t npoints=iPoint++; | |
223 | Int_t index=0; | |
224 | Short_t* mask=new Short_t[npoints]; | |
225 | for(Int_t ip=0;ip<npoints;ip++) mask[ip]=-1; | |
226 | ||
227 | for(Int_t ip1=0;ip1<npoints;ip1++){ | |
228 | if(mask[ip1]==-1) mask[ip1]=index++; | |
229 | for(Int_t ip2=ip1+1; ip2<npoints; ip2++){ | |
230 | if(mask[ip2]==mask[ip1] && mask[ip2]!=-1) continue; | |
231 | Double_t dist2=(xP[ip1]-xP[ip2])*(xP[ip1]-xP[ip2]); | |
232 | dist2+=(yP[ip1]-yP[ip2])*(yP[ip1]-yP[ip2]); | |
233 | dist2+=(zP[ip1]-zP[ip2])*(zP[ip1]-zP[ip2]); | |
234 | if(dist2<fCutOnPairs*fCutOnPairs){ | |
235 | if(mask[ip2]==-1) mask[ip2]=mask[ip1]; | |
236 | else{ | |
237 | for(Int_t ip=0; ip<npoints;ip++){ | |
238 | if(mask[ip]==mask[ip2]) mask[ip]=mask[ip1]; | |
239 | } | |
240 | } | |
241 | } | |
242 | } | |
243 | } | |
244 | ||
245 | ||
246 | // Count multiplicity of trackelts in clusters | |
247 | UInt_t* isIndUsed=new UInt_t[index+1]; | |
248 | for(Int_t ind=0;ind<index+1;ind++) isIndUsed[ind]=0; | |
249 | for(Int_t ip=0; ip<npoints;ip++){ | |
250 | Int_t ind=mask[ip]; | |
251 | isIndUsed[ind]++; | |
252 | } | |
253 | ||
254 | // Count clusters/vertices and sort according to multiplicity | |
255 | Int_t nClusters=0; | |
256 | Int_t* sortedIndex=new Int_t[index+1]; | |
257 | for(Int_t ind1=0;ind1<index+1;ind1++){ | |
258 | if(isIndUsed[ind1]<=1) isIndUsed[ind1]=0; | |
259 | else nClusters++; | |
260 | UInt_t cap=9999999; | |
261 | if(ind1>0) cap=isIndUsed[sortedIndex[ind1-1]]; | |
262 | UInt_t bigger=0; | |
263 | Int_t biggerindex=-1; | |
264 | for(Int_t ind2=0;ind2<index+1;ind2++){ | |
265 | Bool_t use=kTRUE; | |
266 | for(Int_t ind3=0; ind3<ind1; ind3++) | |
267 | if(ind2==sortedIndex[ind3]) use=kFALSE; | |
268 | if(use && isIndUsed[ind2]>bigger && isIndUsed[ind2]<=cap){ | |
269 | bigger=isIndUsed[ind2]; | |
270 | biggerindex=ind2; | |
271 | } | |
272 | } | |
273 | sortedIndex[ind1]=biggerindex; | |
274 | } | |
275 | AliDebug(3,Form("Number of clusters before merging = %d\n",nClusters)); | |
276 | ||
277 | // Assign lines to clusters/vertices and merge clusters which share 1 line | |
278 | Int_t nClustersAfterMerge=nClusters; | |
279 | Int_t* belongsTo=new Int_t[nLines]; | |
280 | for(Int_t ilin=0; ilin<nLines; ilin++) belongsTo[ilin]=-1; | |
281 | for(Int_t iclu=0;iclu<nClusters;iclu++){ | |
282 | Int_t actualCluIndex=iclu; | |
283 | for(Int_t ip=0; ip<npoints;ip++){ | |
284 | if(mask[ip]==sortedIndex[iclu]){ | |
285 | Int_t ind1=index1[ip]; | |
286 | if(belongsTo[ind1]==-1) belongsTo[ind1]=actualCluIndex; | |
287 | else if(belongsTo[ind1]<actualCluIndex){ | |
288 | Int_t newCluIndex=belongsTo[ind1]; | |
289 | for(Int_t ilin=0; ilin<nLines; ilin++){ | |
290 | if(belongsTo[ilin]==actualCluIndex) belongsTo[ilin]=newCluIndex; | |
291 | } | |
292 | AliDebug(10,Form("Merged cluster %d with %d\n",actualCluIndex,newCluIndex)); | |
293 | actualCluIndex=newCluIndex; | |
294 | nClustersAfterMerge--; | |
295 | } | |
296 | Int_t ind2=index2[ip]; | |
297 | if(belongsTo[ind2]==-1) belongsTo[ind2]=actualCluIndex; | |
298 | else if(belongsTo[ind2]<actualCluIndex){ | |
299 | Int_t newCluIndex=belongsTo[ind2]; | |
300 | for(Int_t ilin=0; ilin<nLines; ilin++){ | |
301 | if(belongsTo[ilin]==actualCluIndex) belongsTo[ilin]=newCluIndex; | |
302 | } | |
303 | AliDebug(10,Form("Merged cluster %d with %d\n",actualCluIndex,newCluIndex)); | |
304 | actualCluIndex=newCluIndex; | |
305 | nClustersAfterMerge--; | |
306 | } | |
307 | } | |
308 | } | |
309 | } | |
310 | AliDebug(3,Form("Number of clusters after merging = %d\n",nClustersAfterMerge)); | |
311 | ||
312 | // Count lines associated to each cluster/vertex | |
313 | UInt_t *cluSize=new UInt_t[nClusters]; | |
314 | for(Int_t iclu=0;iclu<nClusters;iclu++){ | |
315 | cluSize[iclu]=0; | |
316 | for(Int_t ilin=0; ilin<nLines; ilin++){ | |
317 | if(belongsTo[ilin]==iclu) cluSize[iclu]++; | |
318 | } | |
319 | } | |
320 | ||
321 | // Count good vertices (>1 associated tracklet) | |
322 | UInt_t nGoodVert=0; | |
323 | for(Int_t iclu=0;iclu<nClusters;iclu++){ | |
324 | AliDebug(3,Form("Vertex %d Size=%d\n",iclu,cluSize[iclu])); | |
325 | if(cluSize[iclu]>1) nGoodVert++; | |
326 | } | |
327 | ||
328 | AliDebug(1,Form("Number of good vertices = %d\n",nGoodVert)); | |
329 | // Calculate vertex coordinates for each cluster | |
330 | if(nGoodVert>0){ | |
331 | fVertArray = new AliESDVertex[nGoodVert]; | |
332 | Int_t iVert=0; | |
333 | for(Int_t iclu=0;iclu<nClusters;iclu++){ | |
334 | Int_t size=cluSize[iclu]; | |
335 | if(size>1){ | |
336 | AliStrLine **arrlin = new AliStrLine*[size]; | |
337 | Int_t nFilled=0; | |
338 | for(Int_t ilin=0; ilin<nLines; ilin++){ | |
339 | if(belongsTo[ilin]==iclu){ | |
340 | arrlin[nFilled++] = dynamic_cast<AliStrLine*>(fLines[ilin]); | |
341 | } | |
342 | } | |
343 | AliDebug(3,Form("Vertex %d N associated tracklets = %d out of %d\n",iVert,size,nFilled)); | |
344 | ||
345 | fVertArray[iVert]=AliVertexerTracks::TrackletVertexFinder(arrlin,nFilled); | |
346 | Double_t peak[3]; | |
347 | fVertArray[iVert].GetXYZ(peak); | |
348 | AliStrLine **arrlin2 = new AliStrLine*[size]; | |
349 | Int_t nFilled2=0; | |
350 | for(Int_t i=0; i<nFilled;i++){ | |
351 | AliStrLine *l1 = arrlin[i]; | |
352 | if(l1->GetDistFromPoint(peak)< fDCAcut) | |
353 | arrlin2[nFilled2++] = dynamic_cast<AliStrLine*>(l1); | |
354 | } | |
355 | if(nFilled2>1){ | |
356 | AliDebug(3,Form("Vertex %d recalculated with %d tracklets\n",iVert,nFilled2)); | |
357 | fVertArray[iVert]=AliVertexerTracks::TrackletVertexFinder(arrlin2,nFilled2); | |
358 | } | |
359 | delete [] arrlin; | |
360 | delete [] arrlin2; | |
361 | ++iVert; | |
362 | } | |
363 | } | |
364 | ||
365 | if(nGoodVert > 1){ | |
366 | fIsPileup = kTRUE; | |
367 | fNTrpuv = fVertArray[1].GetNContributors(); | |
368 | fZpuv = fVertArray[1].GetZv(); | |
369 | } | |
370 | ||
371 | Double_t vRadius=TMath::Sqrt(fVertArray[0].GetXv()*fVertArray[0].GetXv()+fVertArray[0].GetYv()*fVertArray[0].GetYv()); | |
372 | if(vRadius<GetPipeRadius() && fVertArray[0].GetNContributors()>0){ | |
373 | Double_t position[3]={fVertArray[0].GetXv(),fVertArray[0].GetYv(),fVertArray[0].GetZv()}; | |
374 | Double_t covmatrix[6]; | |
375 | fVertArray[0].GetCovMatrix(covmatrix); | |
376 | Double_t chi2=99999.; | |
377 | Int_t nContr=fVertArray[0].GetNContributors(); | |
378 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
379 | fCurrentVertex->SetTitle("vertexer: 3D"); | |
380 | fCurrentVertex->SetName("SPDVertex3D"); | |
381 | fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion()); | |
382 | } | |
383 | } | |
384 | ||
385 | delete [] index1; | |
386 | delete [] index2; | |
387 | delete [] mask; | |
388 | delete [] isIndUsed; | |
389 | delete [] sortedIndex; | |
390 | delete [] belongsTo; | |
391 | delete [] cluSize; | |
392 | delete [] xP; | |
393 | delete [] yP; | |
394 | delete [] zP; | |
395 | } | |
396 | //______________________________________________________________________ | |
397 | void AliITSVertexer3D::FindVertex3DIterativeMM(){ | |
6b4d9537 | 398 | // Defines the AliESDVertex for the current event |
399 | Int_t numsor=fLines.GetEntriesFast()*(fLines.GetEntriesFast()-1)/2; | |
400 | //cout<<"AliITSVertexer3D::FindVertexForCurentEvent: Number of tracklets selected for vertexing "<<fLines.GetEntriesFast()<<"; Number of pairs: "<<numsor<<endl; | |
401 | AliITSSortTrkl srt(fLines,numsor,fCutOnPairs,fCoarseMaxRCut); | |
402 | srt.FindClusters(); | |
403 | AliInfo(Form("Number of vertices: %d",srt.GetNumberOfClusters())); | |
404 | ||
405 | fNoVertices = srt.GetNumberOfClusters(); | |
406 | //printf("fNoVertices = %d \n",fNoVertices); | |
407 | if(fNoVertices>0){ | |
408 | fVertArray = new AliESDVertex[fNoVertices]; | |
409 | for(Int_t kk=0; kk<srt.GetNumberOfClusters(); kk++){ | |
410 | Int_t size = 0; | |
411 | Int_t *labels = srt.GetTrackletsLab(kk,size); | |
412 | /* | |
413 | Int_t *pairs = srt.GetClusters(kk); | |
414 | Int_t nopai = srt.GetSizeOfCluster(kk); | |
415 | cout<<"***** Vertex number "<<kk<<". Pairs: \n"; | |
416 | for(Int_t jj=0;jj<nopai;jj++){ | |
417 | cout<<pairs[jj]<<" - "; | |
418 | if(jj>0 & jj%8==0)cout<<endl; | |
419 | } | |
420 | cout<<endl; | |
421 | cout<<"***** Vertex number "<<kk<<". Labels: \n"; | |
422 | */ | |
423 | AliStrLine **tclo = new AliStrLine* [size]; | |
424 | for(Int_t jj=0;jj<size;jj++){ | |
425 | // cout<<labels[jj]<<" - "; | |
426 | // if(jj>0 & jj%8==0)cout<<endl; | |
427 | tclo[jj] = dynamic_cast<AliStrLine*>(fLines[labels[jj]]); | |
428 | } | |
429 | // cout<<endl; | |
430 | delete []labels; | |
431 | fVertArray[kk]=AliVertexerTracks::TrackletVertexFinder(tclo,size); | |
432 | delete [] tclo; | |
433 | // fVertArray[kk].PrintStatus(); | |
434 | if(kk == 1){ | |
435 | // at least one second vertex is present | |
436 | fIsPileup = kTRUE; | |
437 | fNTrpuv = fVertArray[kk].GetNContributors(); | |
438 | fZpuv = fVertArray[kk].GetZv(); | |
439 | } | |
440 | } | |
0496a544 | 441 | Double_t vRadius=TMath::Sqrt(fVertArray[0].GetXv()*fVertArray[0].GetXv()+fVertArray[0].GetYv()*fVertArray[0].GetYv()); |
6b4d9537 | 442 | if(vRadius<GetPipeRadius() && fVertArray[0].GetNContributors()>0){ |
443 | Double_t position[3]={fVertArray[0].GetXv(),fVertArray[0].GetYv(),fVertArray[0].GetZv()}; | |
444 | Double_t covmatrix[6]; | |
445 | fVertArray[0].GetCovMatrix(covmatrix); | |
446 | Double_t chi2=99999.; | |
447 | Int_t nContr=fVertArray[0].GetNContributors(); | |
448 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
449 | fCurrentVertex->SetTitle("vertexer: 3D"); | |
450 | fCurrentVertex->SetName("SPDVertex3D"); | |
451 | fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion()); | |
452 | } | |
453 | } | |
454 | ||
455 | } | |
456 | ||
457 | //______________________________________________________________________ | |
458 | Bool_t AliITSVertexer3D::DistBetweenVertices(AliESDVertex &a, AliESDVertex &b, Double_t test, Double_t &dist){ | |
459 | // method to compare the distance between vertices a and b with "test" | |
460 | //it returns kTRUE is the distance is less or equal to test | |
461 | dist = (a.GetX()-b.GetX()) * (a.GetX()-b.GetX()); | |
462 | dist += (a.GetY()-b.GetY()) * (a.GetY()-b.GetY()); | |
463 | dist += (a.GetZ()-b.GetZ()) * (a.GetZ()-b.GetZ()); | |
464 | dist = TMath::Sqrt(dist); | |
465 | if(dist <= test)return kTRUE; | |
466 | return kFALSE; | |
467 | } | |
468 | ||
469 | ||
70c95f95 | 470 | //______________________________________________________________________ |
308c2f7c | 471 | Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){ |
70c95f95 | 472 | // All the possible combinations between recpoints on layer 1and 2 are |
473 | // considered. Straight lines (=tracklets)are formed. | |
474 | // The tracklets are processed in Prepare3DVertex | |
70c95f95 | 475 | |
70c95f95 | 476 | TClonesArray *itsRec = 0; |
cd1c2af1 | 477 | if(optCuts==0) fZHisto->Reset(); |
6b4d9537 | 478 | // gc1 are local and global coordinates for layer 1 |
0496a544 | 479 | Float_t gc1f[3]={0.,0.,0.}; |
480 | Double_t gc1[3]={0.,0.,0.}; | |
6b4d9537 | 481 | // gc2 are local and global coordinates for layer 2 |
0496a544 | 482 | Float_t gc2f[3]={0.,0.,0.}; |
483 | Double_t gc2[3]={0.,0.,0.}; | |
b21c1af0 | 484 | AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance(); |
485 | itsRec=rpcont->FetchClusters(0,itsClusterTree); | |
486 | if(!rpcont->IsSPDActive()){ | |
487 | AliWarning("No SPD rec points found, 3D vertex not calculated"); | |
6b4d9537 | 488 | return -1; |
489 | } | |
70c95f95 | 490 | |
491 | // Set values for cuts | |
0496a544 | 492 | Double_t xbeam=GetNominalPos()[0]; |
493 | Double_t ybeam=GetNominalPos()[1]; | |
494 | Double_t zvert=0.; | |
495 | Double_t deltaPhi=fCoarseDiffPhiCut; | |
496 | Double_t deltaR=fCoarseMaxRCut; | |
497 | Double_t dZmax=fZCutDiamond; | |
87104b06 | 498 | if(optCuts==1){ |
27167524 | 499 | xbeam=fVert3D.GetXv(); |
500 | ybeam=fVert3D.GetYv(); | |
501 | zvert=fVert3D.GetZv(); | |
70c95f95 | 502 | deltaPhi = fDiffPhiMax; |
05d1294c | 503 | deltaR=fMaxRCut; |
504 | dZmax=fMaxZCut; | |
87104b06 | 505 | if(fPileupAlgo == 2){ |
506 | dZmax=fZCutDiamond; | |
507 | deltaR=fMaxRCut2; | |
508 | } | |
6b4d9537 | 509 | } else if(optCuts==2){ |
cd1c2af1 | 510 | xbeam=fVert3D.GetXv(); |
511 | ybeam=fVert3D.GetYv(); | |
c238ba33 | 512 | deltaPhi = fDiffPhiforPileup; |
cd1c2af1 | 513 | deltaR=fMaxRCut; |
70c95f95 | 514 | } |
cd1c2af1 | 515 | |
70c95f95 | 516 | Int_t nrpL1 = 0; // number of rec points on layer 1 |
517 | Int_t nrpL2 = 0; // number of rec points on layer 2 | |
1cc75a0b | 518 | nrpL1=rpcont->GetNClustersInLayerFast(1); |
519 | nrpL2=rpcont->GetNClustersInLayerFast(2); | |
70c95f95 | 520 | if(nrpL1 == 0 || nrpL2 == 0){ |
fea0ad94 | 521 | AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",nrpL1,nrpL2)); |
70c95f95 | 522 | return -1; |
523 | } | |
cd706e57 | 524 | AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2)); |
9801e044 | 525 | if(nrpL1>fMaxNumOfCl || nrpL2>fMaxNumOfCl){ |
526 | AliWarning(Form("Too many recpoints on SPD(%d %d ), call vertexerZ",nrpL1,nrpL2)); | |
527 | return -1; | |
528 | } | |
70c95f95 | 529 | |
530 | Double_t a[3]={xbeam,ybeam,0.}; | |
531 | Double_t b[3]={xbeam,ybeam,10.}; | |
532 | AliStrLine zeta(a,b,kTRUE); | |
322e007a | 533 | static Double_t bField=TMath::Abs(AliTracker::GetBz()/10.); //T |
6a0d56b8 | 534 | SetMeanPPtSelTracks(bField); |
87104b06 | 535 | |
70c95f95 | 536 | Int_t nolines = 0; |
537 | // Loop on modules of layer 1 | |
1cc75a0b | 538 | Int_t firstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1); |
539 | Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1; | |
6a0d56b8 | 540 | for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1 |
8c42830a | 541 | if(!fUseModule[modul1]) continue; |
410cfb90 | 542 | UShort_t ladder=modul1/4+1; // ladders are numbered starting from 1 |
b21c1af0 | 543 | TClonesArray *prpl1=rpcont->UncheckedGetClusters(modul1); |
544 | Int_t nrecp1 = prpl1->GetEntries(); | |
70c95f95 | 545 | for(Int_t j=0;j<nrecp1;j++){ |
cd1c2af1 | 546 | if(j>kMaxCluPerMod) continue; |
547 | UShort_t idClu1=modul1*kMaxCluPerMod+j; | |
548 | if(fUsedCluster.TestBitNumber(idClu1)) continue; | |
b21c1af0 | 549 | AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1->At(j); |
0496a544 | 550 | recp1->GetGlobalXYZ(gc1f); |
551 | for(Int_t ico=0;ico<3;ico++)gc1[ico]=gc1f[ico]; | |
552 | ||
6a0d56b8 | 553 | Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam); |
70c95f95 | 554 | if(phi1<0)phi1=2*TMath::Pi()+phi1; |
27167524 | 555 | for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){ |
556 | for(Int_t k=0;k<4;k++){ | |
557 | Int_t ladmod=fLadders[ladder-1]+ladl2; | |
32e63e47 | 558 | if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2); |
559 | Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1); | |
8c42830a | 560 | if(!fUseModule[modul2]) continue; |
b21c1af0 | 561 | itsRec=rpcont->UncheckedGetClusters(modul2); |
27167524 | 562 | Int_t nrecp2 = itsRec->GetEntries(); |
563 | for(Int_t j2=0;j2<nrecp2;j2++){ | |
cd1c2af1 | 564 | if(j2>kMaxCluPerMod) continue; |
565 | UShort_t idClu2=modul2*kMaxCluPerMod+j2; | |
566 | if(fUsedCluster.TestBitNumber(idClu2)) continue; | |
6a0d56b8 | 567 | AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2); |
0496a544 | 568 | recp2->GetGlobalXYZ(gc2f); |
569 | for(Int_t ico=0;ico<3;ico++)gc2[ico]=gc2f[ico]; | |
27167524 | 570 | Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam); |
571 | if(phi2<0)phi2=2*TMath::Pi()+phi2; | |
572 | Double_t diff = TMath::Abs(phi2-phi1); | |
573 | if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; | |
54091147 | 574 | if(optCuts==0 && diff<fDiffPhiforPileup){ |
0496a544 | 575 | Double_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); |
576 | Double_t zc1=gc1[2]; | |
577 | Double_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); | |
578 | Double_t zc2=gc2[2]; | |
579 | Double_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius | |
cd1c2af1 | 580 | fZHisto->Fill(zr0); |
581 | } | |
27167524 | 582 | if(diff>deltaPhi)continue; |
6a0d56b8 | 583 | AliStrLine line(gc1,gc2,kTRUE); |
27167524 | 584 | Double_t cp[3]; |
585 | Int_t retcode = line.Cross(&zeta,cp); | |
586 | if(retcode<0)continue; | |
587 | Double_t dca = line.GetDCA(&zeta); | |
588 | if(dca<0.) continue; | |
589 | if(dca>deltaR)continue; | |
590 | Double_t deltaZ=cp[2]-zvert; | |
591 | if(TMath::Abs(deltaZ)>dZmax)continue; | |
6a0d56b8 | 592 | |
6a0d56b8 | 593 | if(nolines == 0){ |
f5f6da22 | 594 | if(fLines.GetEntriesFast()>0)fLines.Clear("C"); |
6a0d56b8 | 595 | } |
596 | Float_t cov[6]; | |
597 | recp2->GetGlobalCov(cov); | |
598 | ||
cd1c2af1 | 599 | |
0496a544 | 600 | Double_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); |
601 | Double_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); | |
602 | Double_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction | |
6a0d56b8 | 603 | |
0496a544 | 604 | Double_t curvErr=0; |
6a0d56b8 | 605 | if(bField>0.00001){ |
0496a544 | 606 | Double_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm |
440198c1 | 607 | Double_t dRad=TMath::Sqrt((gc1[0]-gc2[0])*(gc1[0]-gc2[0])+(gc1[1]-gc2[1])*(gc1[1]-gc2[1])); |
0496a544 | 608 | Double_t aux=dRad/2.+rad1; |
6a0d56b8 | 609 | curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm |
610 | } | |
0496a544 | 611 | Double_t sigmasq[3]; |
6a0d56b8 | 612 | sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor; |
613 | sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor; | |
614 | sigmasq[2]=cov[5]*factor*factor; | |
615 | ||
616 | // Multiple scattering | |
0496a544 | 617 | Double_t pOverMass=fMeanPSelTrk/0.140; |
618 | Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass); | |
619 | Double_t p2=fMeanPSelTrk*fMeanPSelTrk; | |
620 | Double_t rBP=GetPipeRadius(); | |
621 | Double_t dBP=0.08/35.3; // 800 um of Be | |
622 | Double_t dL1=0.01; //approx. 1% of radiation length | |
623 | Double_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP; | |
624 | Double_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1; | |
625 | Double_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1)); | |
626 | Double_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP)); | |
6a0d56b8 | 627 | for(Int_t ico=0; ico<3;ico++){ |
5951029f | 628 | sigmasq[ico]+=rtantheta1*rtantheta1*factor*factor/3.; |
629 | sigmasq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.; | |
6a0d56b8 | 630 | } |
0496a544 | 631 | Double_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.}; |
6a0d56b8 | 632 | if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0]; |
633 | if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1]; | |
634 | if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2]; | |
cd1c2af1 | 635 | new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE,idClu1,idClu2); |
6a0d56b8 | 636 | |
27167524 | 637 | } |
70c95f95 | 638 | } |
70c95f95 | 639 | } |
640 | } | |
70c95f95 | 641 | } |
642 | if(nolines == 0)return -2; | |
643 | return nolines; | |
644 | } | |
645 | ||
70c95f95 | 646 | //______________________________________________________________________ |
05d1294c | 647 | Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ |
70c95f95 | 648 | // Finds the 3D vertex information using tracklets |
649 | Int_t retcode = -1; | |
650 | ||
0496a544 | 651 | Double_t xbeam=GetNominalPos()[0]; |
652 | Double_t ybeam=GetNominalPos()[1]; | |
653 | Double_t zvert=0.; | |
654 | Double_t deltaR=fCoarseMaxRCut; | |
0496a544 | 655 | Double_t dZmax=fZCutDiamond; |
cd1c2af1 | 656 | if(optCuts==1){ |
27167524 | 657 | xbeam=fVert3D.GetXv(); |
658 | ybeam=fVert3D.GetYv(); | |
659 | zvert=fVert3D.GetZv(); | |
05d1294c | 660 | deltaR=fMaxRCut; |
661 | dZmax=fMaxZCut; | |
87104b06 | 662 | if(fPileupAlgo == 2){ |
663 | dZmax=fZCutDiamond; | |
664 | deltaR=fMaxRCut2; | |
665 | } | |
cd1c2af1 | 666 | }else if(optCuts==2){ |
667 | xbeam=fVert3D.GetXv(); | |
668 | ybeam=fVert3D.GetYv(); | |
669 | deltaR=fMaxRCut; | |
05d1294c | 670 | } |
671 | ||
0496a544 | 672 | Double_t rl=-fCoarseMaxRCut; |
673 | Double_t rh=fCoarseMaxRCut; | |
674 | Double_t zl=-fZCutDiamond; | |
675 | Double_t zh=fZCutDiamond; | |
7d766abb | 676 | Int_t nbr=(Int_t)((rh-rl)/fBinSizeR+0.0001); |
677 | Int_t nbz=(Int_t)((zh-zl)/fBinSizeZ+0.0001); | |
678 | Int_t nbrcs=(Int_t)((rh-rl)/(fBinSizeR*2.)+0.0001); | |
679 | Int_t nbzcs=(Int_t)((zh-zl)/(fBinSizeZ*2.)+0.0001); | |
680 | ||
87104b06 | 681 | TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh); |
682 | TH3F *h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh); | |
683 | ||
70c95f95 | 684 | // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed) |
f5f6da22 | 685 | Int_t *validate = new Int_t [fLines.GetEntriesFast()]; |
686 | for(Int_t i=0; i<fLines.GetEntriesFast();i++)validate[i]=0; | |
687 | for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){ | |
f5f6da22 | 688 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); |
689 | for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){ | |
690 | AliStrLine *l2 = (AliStrLine*)fLines.At(j); | |
0496a544 | 691 | Double_t dca=l1->GetDCA(l2); |
05d1294c | 692 | if(dca > fDCAcut || dca<0.00001) continue; |
70c95f95 | 693 | Double_t point[3]; |
694 | Int_t retc = l1->Cross(l2,point); | |
695 | if(retc<0)continue; | |
6a0d56b8 | 696 | Double_t deltaZ=point[2]-zvert; |
697 | if(TMath::Abs(deltaZ)>dZmax)continue; | |
70c95f95 | 698 | Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]); |
05d1294c | 699 | if(rad>fCoarseMaxRCut)continue; |
700 | Double_t deltaX=point[0]-xbeam; | |
701 | Double_t deltaY=point[1]-ybeam; | |
05d1294c | 702 | Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY); |
05d1294c | 703 | if(raddist>deltaR)continue; |
27167524 | 704 | validate[i]=1; |
705 | validate[j]=1; | |
87104b06 | 706 | h3d->Fill(point[0],point[1],point[2]); |
707 | h3dcs->Fill(point[0],point[1],point[2]); | |
70c95f95 | 708 | } |
709 | } | |
710 | ||
711 | ||
712 | ||
70c95f95 | 713 | Int_t numbtracklets=0; |
f5f6da22 | 714 | for(Int_t i=0; i<fLines.GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++; |
6b4d9537 | 715 | if(numbtracklets<2){ |
716 | delete [] validate; | |
87104b06 | 717 | delete h3d; |
718 | delete h3dcs; | |
6b4d9537 | 719 | return retcode; |
720 | } | |
70c95f95 | 721 | |
f5f6da22 | 722 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
723 | if(validate[i]<1)fLines.RemoveAt(i); | |
70c95f95 | 724 | } |
f5f6da22 | 725 | fLines.Compress(); |
726 | AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast())); | |
70c95f95 | 727 | delete [] validate; |
728 | ||
87104b06 | 729 | // Exit here if Pileup Algorithm 2 has been chosen during second loop |
730 | if(fPileupAlgo == 2 && optCuts==1){ | |
731 | delete h3d; | |
732 | delete h3dcs; | |
733 | return 0; | |
734 | } | |
6b4d9537 | 735 | |
7340864d | 736 | // Find peaks in histos |
05d1294c | 737 | |
70c95f95 | 738 | Double_t peak[3]={0.,0.,0.}; |
7340864d | 739 | Int_t ntrkl,ntimes; |
740 | FindPeaks(h3d,peak,ntrkl,ntimes); | |
70c95f95 | 741 | delete h3d; |
0496a544 | 742 | Double_t binsizer=(rh-rl)/nbr; |
743 | Double_t binsizez=(zh-zl)/nbz; | |
7d766abb | 744 | if(optCuts==0 && (ntrkl<=2 || ntimes>1)){ |
7340864d | 745 | ntrkl=0; |
746 | ntimes=0; | |
747 | FindPeaks(h3dcs,peak,ntrkl,ntimes); | |
748 | binsizer=(rh-rl)/nbrcs; | |
749 | binsizez=(zh-zl)/nbzcs; | |
750 | if(ntrkl==1 || ntimes>1){delete h3dcs; return retcode;} | |
751 | } | |
752 | delete h3dcs; | |
753 | ||
70c95f95 | 754 | // Second selection loop |
7340864d | 755 | |
0496a544 | 756 | Double_t bs=(binsizer+binsizez)/2.; |
f5f6da22 | 757 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
758 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); | |
759 | if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i); | |
70c95f95 | 760 | } |
f5f6da22 | 761 | fLines.Compress(); |
762 | AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast())); | |
70c95f95 | 763 | |
f5f6da22 | 764 | if(fLines.GetEntriesFast()>1){ |
2ecf7a6f | 765 | retcode=0; |
27167524 | 766 | // find a first candidate for the primary vertex |
f5f6da22 | 767 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); |
70c95f95 | 768 | // make a further selection on tracklets based on this first candidate |
27167524 | 769 | fVert3D.GetXYZ(peak); |
770 | AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2])); | |
410cfb90 | 771 | Int_t *validate2 = new Int_t [fLines.GetEntriesFast()]; |
772 | for(Int_t i=0; i<fLines.GetEntriesFast();i++) validate2[i]=1; | |
f5f6da22 | 773 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
410cfb90 | 774 | if(validate2[i]==0) continue; |
f5f6da22 | 775 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); |
776 | if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i); | |
410cfb90 | 777 | if(optCuts==2){ // temporarily only for pileup |
778 | for(Int_t j=i+1; j<fLines.GetEntriesFast();j++){ | |
779 | AliStrLine *l2 = (AliStrLine*)fLines.At(j); | |
780 | if(l1->GetDCA(l2)<0.00001){ | |
781 | Int_t delta1=(Int_t)l1->GetIdPoint(0)-(Int_t)l2->GetIdPoint(0); | |
782 | Int_t delta2=(Int_t)l1->GetIdPoint(1)-(Int_t)l2->GetIdPoint(1); | |
783 | Int_t deltamod1=(Int_t)l1->GetIdPoint(0)/kMaxCluPerMod | |
784 | -(Int_t)l2->GetIdPoint(0)/kMaxCluPerMod; | |
785 | Int_t deltamod2=(Int_t)l1->GetIdPoint(1)/kMaxCluPerMod | |
786 | -(Int_t)l2->GetIdPoint(1)/kMaxCluPerMod; | |
787 | // remove tracklets sharing a point | |
788 | if( (delta1==0 && deltamod2==0) || | |
789 | (delta2==0 && deltamod1==0) ) validate2[j]=0; | |
790 | } | |
791 | } | |
792 | } | |
70c95f95 | 793 | } |
410cfb90 | 794 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
795 | if(validate2[i]==0) fLines.RemoveAt(i); | |
796 | } | |
797 | delete [] validate2; | |
f5f6da22 | 798 | fLines.Compress(); |
799 | AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast())); | |
2ecf7a6f | 800 | if(fLines.GetEntriesFast()>1){// this new tracklet selection is used |
801 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); | |
802 | } | |
70c95f95 | 803 | } |
804 | return retcode; | |
805 | } | |
806 | ||
6a0d56b8 | 807 | //________________________________________________________ |
0496a544 | 808 | void AliITSVertexer3D::SetMeanPPtSelTracks(Double_t fieldTesla){ |
5951029f | 809 | // Sets mean values of Pt based on the field |
810 | // for P (used in multiple scattering) the most probable value is used | |
6a0d56b8 | 811 | if(TMath::Abs(fieldTesla-0.5)<0.01){ |
5951029f | 812 | SetMeanPSelTracks(0.375); |
6a0d56b8 | 813 | SetMeanPtSelTracks(0.630); |
814 | }else if(TMath::Abs(fieldTesla-0.4)<0.01){ | |
5951029f | 815 | SetMeanPSelTracks(0.375); |
6a0d56b8 | 816 | SetMeanPtSelTracks(0.580); |
817 | }else if(TMath::Abs(fieldTesla-0.2)<0.01){ | |
5951029f | 818 | SetMeanPSelTracks(0.375); |
6a0d56b8 | 819 | SetMeanPtSelTracks(0.530); |
820 | }else if(fieldTesla<0.00001){ | |
5951029f | 821 | SetMeanPSelTracks(0.375); |
822 | SetMeanPtSelTracks(0.230); | |
6a0d56b8 | 823 | }else{ |
824 | SetMeanPSelTracks(); | |
825 | SetMeanPtSelTracks(); | |
70c95f95 | 826 | } |
70c95f95 | 827 | } |
828 | ||
7340864d | 829 | //________________________________________________________ |
830 | void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){ | |
831 | // Finds bin with max contents in 3D histo of tracket intersections | |
832 | TAxis *xax = histo->GetXaxis(); | |
833 | TAxis *yax = histo->GetYaxis(); | |
834 | TAxis *zax = histo->GetZaxis(); | |
835 | peak[0]=0.; | |
836 | peak[1]=0.; | |
837 | peak[2]=0.; | |
838 | nOfTracklets = 0; | |
839 | nOfTimes=0; | |
7d766abb | 840 | Int_t peakbin[3]={0,0,0}; |
841 | Int_t peak2bin[3]={-1,-1,-1}; | |
842 | Int_t bc2=-1; | |
7340864d | 843 | for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){ |
242c8025 | 844 | Double_t xval = xax->GetBinCenter(i); |
7340864d | 845 | for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){ |
242c8025 | 846 | Double_t yval = yax->GetBinCenter(j); |
7340864d | 847 | for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){ |
242c8025 | 848 | Double_t zval = zax->GetBinCenter(k); |
7340864d | 849 | Int_t bc =(Int_t)histo->GetBinContent(i,j,k); |
7d766abb | 850 | if(bc==0) continue; |
7340864d | 851 | if(bc>nOfTracklets){ |
7d766abb | 852 | nOfTracklets=bc; |
7340864d | 853 | peak[2] = zval; |
854 | peak[1] = yval; | |
855 | peak[0] = xval; | |
7d766abb | 856 | peakbin[2] = k; |
857 | peakbin[1] = j; | |
858 | peakbin[0] = i; | |
859 | peak2bin[2] = -1; | |
860 | peak2bin[1] = -1; | |
861 | peak2bin[0] = -1; | |
862 | bc2=-1; | |
7340864d | 863 | nOfTimes = 1; |
06a7cbee | 864 | }else if(bc==nOfTracklets){ |
7d766abb | 865 | if(TMath::Abs(i-peakbin[0])<=1 && TMath::Abs(j-peakbin[1])<=1 && TMath::Abs(k-peakbin[2])<=1){ |
866 | peak2bin[2] = k; | |
867 | peak2bin[1] = j; | |
868 | peak2bin[0] = i; | |
869 | bc2=bc; | |
870 | nOfTimes = 1; | |
871 | }else{ | |
872 | nOfTimes++; | |
873 | } | |
7340864d | 874 | } |
875 | } | |
876 | } | |
877 | } | |
7d766abb | 878 | if(peak2bin[0]>=-1 && bc2!=-1){ // two contiguous peak-cells with same contents |
879 | peak[0]=0.5*(xax->GetBinCenter(peakbin[0])+xax->GetBinCenter(peak2bin[0])); | |
880 | peak[1]=0.5*(yax->GetBinCenter(peakbin[1])+yax->GetBinCenter(peak2bin[1])); | |
881 | peak[2]=0.5*(zax->GetBinCenter(peakbin[2])+zax->GetBinCenter(peak2bin[2])); | |
882 | nOfTracklets+=bc2; | |
883 | nOfTimes=1; | |
884 | } | |
7340864d | 885 | } |
cd1c2af1 | 886 | //________________________________________________________ |
887 | void AliITSVertexer3D::MarkUsedClusters(){ | |
888 | // Mark clusters of tracklets used in vertex claulation | |
889 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ | |
890 | AliStrLine *lin = (AliStrLine*)fLines.At(i); | |
891 | Int_t idClu1=lin->GetIdPoint(0); | |
892 | Int_t idClu2=lin->GetIdPoint(1); | |
893 | fUsedCluster.SetBitNumber(idClu1); | |
894 | fUsedCluster.SetBitNumber(idClu2); | |
895 | } | |
896 | } | |
897 | //________________________________________________________ | |
898 | Int_t AliITSVertexer3D::RemoveTracklets(){ | |
899 | // Remove trackelts close to first found vertex | |
900 | Double_t vert[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()}; | |
901 | Int_t nRemoved=0; | |
902 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ | |
903 | AliStrLine *lin = (AliStrLine*)fLines.At(i); | |
904 | if(lin->GetDistFromPoint(vert)<fDCAforPileup){ | |
905 | Int_t idClu1=lin->GetIdPoint(0); | |
906 | Int_t idClu2=lin->GetIdPoint(1); | |
907 | fUsedCluster.SetBitNumber(idClu1); | |
908 | fUsedCluster.SetBitNumber(idClu2); | |
909 | fLines.RemoveAt(i); | |
910 | ++nRemoved; | |
911 | } | |
912 | } | |
913 | fLines.Compress(); | |
914 | return nRemoved; | |
915 | } | |
916 | //________________________________________________________ | |
917 | void AliITSVertexer3D::FindOther3DVertices(TTree *itsClusterTree){ | |
918 | // pileup identification based on 3D vertexing with not used clusters | |
919 | MarkUsedClusters(); | |
920 | fLines.Clear("C"); | |
921 | Int_t nolines = FindTracklets(itsClusterTree,2); | |
922 | if(nolines>=2){ | |
923 | Int_t nr=RemoveTracklets(); | |
924 | nolines-=nr; | |
925 | if(nolines>=2){ | |
926 | Int_t rc=Prepare3DVertex(2); | |
927 | if(rc==0){ | |
928 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); | |
929 | if(fVert3D.GetNContributors()>=fMinTrackletsForPilup){ | |
930 | fIsPileup=kTRUE; | |
6b4d9537 | 931 | fNoVertices=2; |
932 | fVertArray = new AliESDVertex[2]; | |
933 | fVertArray[0]=(*fCurrentVertex); | |
934 | fVertArray[1]=fVert3D; | |
cd1c2af1 | 935 | fZpuv=fVert3D.GetZv(); |
936 | fNTrpuv=fVert3D.GetNContributors(); | |
937 | } | |
938 | } | |
939 | } | |
940 | } | |
941 | } | |
942 | //______________________________________________________________________ | |
943 | void AliITSVertexer3D::PileupFromZ(){ | |
944 | // Calls the pileup algorithm of ALiITSVertexerZ | |
945 | Int_t binmin, binmax; | |
946 | Int_t nPeaks=AliITSVertexerZ::GetPeakRegion(fZHisto,binmin,binmax); | |
947 | if(nPeaks==2)AliWarning("2 peaks found"); | |
948 | Int_t firstPeakCont=0; | |
242c8025 | 949 | Double_t firstPeakPos=0.; |
cd1c2af1 | 950 | for(Int_t i=binmin-1;i<=binmax+1;i++){ |
6b4d9537 | 951 | firstPeakCont+=static_cast<Int_t>(fZHisto->GetBinContent(i)); |
cd1c2af1 | 952 | firstPeakPos+=fZHisto->GetBinContent(i)*fZHisto->GetBinCenter(i); |
953 | } | |
73315c4a | 954 | if(firstPeakCont>0){ |
955 | firstPeakPos/=firstPeakCont; | |
956 | Int_t ncontr2=0; | |
957 | if(firstPeakCont>fMinTrackletsForPilup){ | |
958 | Float_t secPeakPos; | |
959 | ncontr2=AliITSVertexerZ::FindSecondPeak(fZHisto,binmin,binmax,secPeakPos); | |
960 | if(ncontr2>=fMinTrackletsForPilup){ | |
961 | fIsPileup=kTRUE; | |
6b4d9537 | 962 | fNoVertices=2; |
963 | AliESDVertex secondVert(secPeakPos,0.1,ncontr2); | |
964 | fVertArray = new AliESDVertex[2]; | |
965 | fVertArray[0]=(*fCurrentVertex); | |
966 | fVertArray[1]=secondVert; | |
73315c4a | 967 | fZpuv=secPeakPos; |
968 | fNTrpuv=ncontr2; | |
969 | } | |
cd1c2af1 | 970 | } |
971 | } | |
972 | } | |
70c95f95 | 973 | //________________________________________________________ |
974 | void AliITSVertexer3D::PrintStatus() const { | |
975 | // Print current status | |
87104b06 | 976 | printf("========= First step selections =====================\n"); |
6a0d56b8 | 977 | printf("Cut on diamond (Z) %f\n",fZCutDiamond); |
87104b06 | 978 | printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut); |
979 | printf("Loose cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut); | |
6a0d56b8 | 980 | printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut); |
87104b06 | 981 | printf("========= Second step selections ====================\n"); |
982 | printf("Cut on tracklet-to-first-vertex Z distance %f\n",fMaxZCut); | |
cd1c2af1 | 983 | printf("Max Phi difference: %f\n",fDiffPhiMax); |
87104b06 | 984 | printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut); |
985 | printf("Cut on tracklet DCA to beam axis (algo2) %f\n",fMaxRCut2); | |
986 | printf("========= Pileup selections =========================\n"); | |
cd1c2af1 | 987 | printf("Pileup algo: %d\n",fPileupAlgo); |
87104b06 | 988 | printf("Min DCA to 1st vertex for pileup (algo 0 and 1): %f\n",fDCAforPileup); |
989 | printf("Cut on distance between pair-vertices (algo 2): %f\n",fCutOnPairs); | |
9801e044 | 990 | printf("Maximum number of clusters allowed on L1 or L2: %d\n",fMaxNumOfCl); |
6a0d56b8 | 991 | printf("=======================================================\n"); |
70c95f95 | 992 | } |
410cfb90 | 993 |