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 |
518 | |
519 | // By default irstL1=0 and lastL1=79 |
6a0d56b8 |
520 | Int_t firstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1); |
32e63e47 |
521 | Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1; |
6a0d56b8 |
522 | for(Int_t module= firstL1; module<=lastL1;module++){ // count number of recopints on layer 1 |
b21c1af0 |
523 | itsRec=rpcont->UncheckedGetClusters(module); |
70c95f95 |
524 | nrpL1+= itsRec->GetEntries(); |
70c95f95 |
525 | } |
6a0d56b8 |
526 | //By default firstL2=80 and lastL2=239 |
527 | Int_t firstL2 = AliITSgeomTGeo::GetModuleIndex(2,1,1); |
32e63e47 |
528 | Int_t lastL2 = AliITSgeomTGeo::GetModuleIndex(3,1,1)-1; |
6a0d56b8 |
529 | for(Int_t module= firstL2; module<=lastL2;module++){ // count number of recopints on layer 2 |
b21c1af0 |
530 | itsRec=rpcont->UncheckedGetClusters(module); |
70c95f95 |
531 | nrpL2+= itsRec->GetEntries(); |
70c95f95 |
532 | } |
533 | if(nrpL1 == 0 || nrpL2 == 0){ |
fea0ad94 |
534 | AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",nrpL1,nrpL2)); |
70c95f95 |
535 | return -1; |
536 | } |
cd706e57 |
537 | AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2)); |
9801e044 |
538 | if(nrpL1>fMaxNumOfCl || nrpL2>fMaxNumOfCl){ |
539 | AliWarning(Form("Too many recpoints on SPD(%d %d ), call vertexerZ",nrpL1,nrpL2)); |
540 | return -1; |
541 | } |
70c95f95 |
542 | |
543 | Double_t a[3]={xbeam,ybeam,0.}; |
544 | Double_t b[3]={xbeam,ybeam,10.}; |
545 | AliStrLine zeta(a,b,kTRUE); |
322e007a |
546 | static Double_t bField=TMath::Abs(AliTracker::GetBz()/10.); //T |
6a0d56b8 |
547 | SetMeanPPtSelTracks(bField); |
87104b06 |
548 | |
70c95f95 |
549 | Int_t nolines = 0; |
550 | // Loop on modules of layer 1 |
6a0d56b8 |
551 | for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1 |
8c42830a |
552 | if(!fUseModule[modul1]) continue; |
27167524 |
553 | UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1 |
b21c1af0 |
554 | TClonesArray *prpl1=rpcont->UncheckedGetClusters(modul1); |
555 | Int_t nrecp1 = prpl1->GetEntries(); |
70c95f95 |
556 | for(Int_t j=0;j<nrecp1;j++){ |
cd1c2af1 |
557 | if(j>kMaxCluPerMod) continue; |
558 | UShort_t idClu1=modul1*kMaxCluPerMod+j; |
559 | if(fUsedCluster.TestBitNumber(idClu1)) continue; |
b21c1af0 |
560 | AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1->At(j); |
0496a544 |
561 | recp1->GetGlobalXYZ(gc1f); |
562 | for(Int_t ico=0;ico<3;ico++)gc1[ico]=gc1f[ico]; |
563 | |
6a0d56b8 |
564 | Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam); |
70c95f95 |
565 | if(phi1<0)phi1=2*TMath::Pi()+phi1; |
27167524 |
566 | for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){ |
567 | for(Int_t k=0;k<4;k++){ |
568 | Int_t ladmod=fLadders[ladder-1]+ladl2; |
32e63e47 |
569 | if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2); |
570 | Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1); |
8c42830a |
571 | if(!fUseModule[modul2]) continue; |
b21c1af0 |
572 | itsRec=rpcont->UncheckedGetClusters(modul2); |
27167524 |
573 | Int_t nrecp2 = itsRec->GetEntries(); |
574 | for(Int_t j2=0;j2<nrecp2;j2++){ |
cd1c2af1 |
575 | if(j2>kMaxCluPerMod) continue; |
576 | UShort_t idClu2=modul2*kMaxCluPerMod+j2; |
577 | if(fUsedCluster.TestBitNumber(idClu2)) continue; |
6a0d56b8 |
578 | AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2); |
0496a544 |
579 | recp2->GetGlobalXYZ(gc2f); |
580 | for(Int_t ico=0;ico<3;ico++)gc2[ico]=gc2f[ico]; |
27167524 |
581 | Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam); |
582 | if(phi2<0)phi2=2*TMath::Pi()+phi2; |
583 | Double_t diff = TMath::Abs(phi2-phi1); |
584 | if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; |
54091147 |
585 | if(optCuts==0 && diff<fDiffPhiforPileup){ |
0496a544 |
586 | Double_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); |
587 | Double_t zc1=gc1[2]; |
588 | Double_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); |
589 | Double_t zc2=gc2[2]; |
590 | Double_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius |
cd1c2af1 |
591 | fZHisto->Fill(zr0); |
592 | } |
27167524 |
593 | if(diff>deltaPhi)continue; |
6a0d56b8 |
594 | AliStrLine line(gc1,gc2,kTRUE); |
27167524 |
595 | Double_t cp[3]; |
596 | Int_t retcode = line.Cross(&zeta,cp); |
597 | if(retcode<0)continue; |
598 | Double_t dca = line.GetDCA(&zeta); |
599 | if(dca<0.) continue; |
600 | if(dca>deltaR)continue; |
601 | Double_t deltaZ=cp[2]-zvert; |
602 | if(TMath::Abs(deltaZ)>dZmax)continue; |
6a0d56b8 |
603 | |
6a0d56b8 |
604 | if(nolines == 0){ |
f5f6da22 |
605 | if(fLines.GetEntriesFast()>0)fLines.Clear("C"); |
6a0d56b8 |
606 | } |
607 | Float_t cov[6]; |
608 | recp2->GetGlobalCov(cov); |
609 | |
cd1c2af1 |
610 | |
0496a544 |
611 | Double_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); |
612 | Double_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); |
613 | Double_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction |
6a0d56b8 |
614 | |
0496a544 |
615 | Double_t curvErr=0; |
6a0d56b8 |
616 | if(bField>0.00001){ |
0496a544 |
617 | Double_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm |
440198c1 |
618 | Double_t dRad=TMath::Sqrt((gc1[0]-gc2[0])*(gc1[0]-gc2[0])+(gc1[1]-gc2[1])*(gc1[1]-gc2[1])); |
0496a544 |
619 | Double_t aux=dRad/2.+rad1; |
6a0d56b8 |
620 | curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm |
621 | } |
0496a544 |
622 | Double_t sigmasq[3]; |
6a0d56b8 |
623 | sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor; |
624 | sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor; |
625 | sigmasq[2]=cov[5]*factor*factor; |
626 | |
627 | // Multiple scattering |
0496a544 |
628 | Double_t pOverMass=fMeanPSelTrk/0.140; |
629 | Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass); |
630 | Double_t p2=fMeanPSelTrk*fMeanPSelTrk; |
631 | Double_t rBP=GetPipeRadius(); |
632 | Double_t dBP=0.08/35.3; // 800 um of Be |
633 | Double_t dL1=0.01; //approx. 1% of radiation length |
634 | Double_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP; |
635 | Double_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1; |
636 | Double_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1)); |
637 | Double_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP)); |
6a0d56b8 |
638 | for(Int_t ico=0; ico<3;ico++){ |
5951029f |
639 | sigmasq[ico]+=rtantheta1*rtantheta1*factor*factor/3.; |
640 | sigmasq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.; |
6a0d56b8 |
641 | } |
0496a544 |
642 | Double_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.}; |
6a0d56b8 |
643 | if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0]; |
644 | if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1]; |
645 | if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2]; |
cd1c2af1 |
646 | new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE,idClu1,idClu2); |
6a0d56b8 |
647 | |
27167524 |
648 | } |
70c95f95 |
649 | } |
70c95f95 |
650 | } |
651 | } |
70c95f95 |
652 | } |
653 | if(nolines == 0)return -2; |
654 | return nolines; |
655 | } |
656 | |
70c95f95 |
657 | //______________________________________________________________________ |
05d1294c |
658 | Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ |
70c95f95 |
659 | // Finds the 3D vertex information using tracklets |
660 | Int_t retcode = -1; |
661 | |
0496a544 |
662 | Double_t xbeam=GetNominalPos()[0]; |
663 | Double_t ybeam=GetNominalPos()[1]; |
664 | Double_t zvert=0.; |
665 | Double_t deltaR=fCoarseMaxRCut; |
0496a544 |
666 | Double_t dZmax=fZCutDiamond; |
cd1c2af1 |
667 | if(optCuts==1){ |
27167524 |
668 | xbeam=fVert3D.GetXv(); |
669 | ybeam=fVert3D.GetYv(); |
670 | zvert=fVert3D.GetZv(); |
05d1294c |
671 | deltaR=fMaxRCut; |
672 | dZmax=fMaxZCut; |
87104b06 |
673 | if(fPileupAlgo == 2){ |
674 | dZmax=fZCutDiamond; |
675 | deltaR=fMaxRCut2; |
676 | } |
cd1c2af1 |
677 | }else if(optCuts==2){ |
678 | xbeam=fVert3D.GetXv(); |
679 | ybeam=fVert3D.GetYv(); |
680 | deltaR=fMaxRCut; |
05d1294c |
681 | } |
682 | |
0496a544 |
683 | Double_t rl=-fCoarseMaxRCut; |
684 | Double_t rh=fCoarseMaxRCut; |
685 | Double_t zl=-fZCutDiamond; |
686 | Double_t zh=fZCutDiamond; |
7d766abb |
687 | Int_t nbr=(Int_t)((rh-rl)/fBinSizeR+0.0001); |
688 | Int_t nbz=(Int_t)((zh-zl)/fBinSizeZ+0.0001); |
689 | Int_t nbrcs=(Int_t)((rh-rl)/(fBinSizeR*2.)+0.0001); |
690 | Int_t nbzcs=(Int_t)((zh-zl)/(fBinSizeZ*2.)+0.0001); |
691 | |
87104b06 |
692 | TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh); |
693 | TH3F *h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh); |
694 | |
70c95f95 |
695 | // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed) |
f5f6da22 |
696 | Int_t *validate = new Int_t [fLines.GetEntriesFast()]; |
697 | for(Int_t i=0; i<fLines.GetEntriesFast();i++)validate[i]=0; |
698 | for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){ |
f5f6da22 |
699 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); |
700 | for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){ |
701 | AliStrLine *l2 = (AliStrLine*)fLines.At(j); |
0496a544 |
702 | Double_t dca=l1->GetDCA(l2); |
05d1294c |
703 | if(dca > fDCAcut || dca<0.00001) continue; |
70c95f95 |
704 | Double_t point[3]; |
705 | Int_t retc = l1->Cross(l2,point); |
706 | if(retc<0)continue; |
6a0d56b8 |
707 | Double_t deltaZ=point[2]-zvert; |
708 | if(TMath::Abs(deltaZ)>dZmax)continue; |
70c95f95 |
709 | Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]); |
05d1294c |
710 | if(rad>fCoarseMaxRCut)continue; |
711 | Double_t deltaX=point[0]-xbeam; |
712 | Double_t deltaY=point[1]-ybeam; |
05d1294c |
713 | Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY); |
05d1294c |
714 | if(raddist>deltaR)continue; |
27167524 |
715 | validate[i]=1; |
716 | validate[j]=1; |
87104b06 |
717 | h3d->Fill(point[0],point[1],point[2]); |
718 | h3dcs->Fill(point[0],point[1],point[2]); |
70c95f95 |
719 | } |
720 | } |
721 | |
722 | |
723 | |
70c95f95 |
724 | Int_t numbtracklets=0; |
f5f6da22 |
725 | for(Int_t i=0; i<fLines.GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++; |
6b4d9537 |
726 | if(numbtracklets<2){ |
727 | delete [] validate; |
87104b06 |
728 | delete h3d; |
729 | delete h3dcs; |
6b4d9537 |
730 | return retcode; |
731 | } |
70c95f95 |
732 | |
f5f6da22 |
733 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
734 | if(validate[i]<1)fLines.RemoveAt(i); |
70c95f95 |
735 | } |
f5f6da22 |
736 | fLines.Compress(); |
737 | AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast())); |
70c95f95 |
738 | delete [] validate; |
739 | |
87104b06 |
740 | // Exit here if Pileup Algorithm 2 has been chosen during second loop |
741 | if(fPileupAlgo == 2 && optCuts==1){ |
742 | delete h3d; |
743 | delete h3dcs; |
744 | return 0; |
745 | } |
6b4d9537 |
746 | |
7340864d |
747 | // Find peaks in histos |
05d1294c |
748 | |
70c95f95 |
749 | Double_t peak[3]={0.,0.,0.}; |
7340864d |
750 | Int_t ntrkl,ntimes; |
751 | FindPeaks(h3d,peak,ntrkl,ntimes); |
70c95f95 |
752 | delete h3d; |
0496a544 |
753 | Double_t binsizer=(rh-rl)/nbr; |
754 | Double_t binsizez=(zh-zl)/nbz; |
7d766abb |
755 | if(optCuts==0 && (ntrkl<=2 || ntimes>1)){ |
7340864d |
756 | ntrkl=0; |
757 | ntimes=0; |
758 | FindPeaks(h3dcs,peak,ntrkl,ntimes); |
759 | binsizer=(rh-rl)/nbrcs; |
760 | binsizez=(zh-zl)/nbzcs; |
761 | if(ntrkl==1 || ntimes>1){delete h3dcs; return retcode;} |
762 | } |
763 | delete h3dcs; |
764 | |
70c95f95 |
765 | // Second selection loop |
7340864d |
766 | |
0496a544 |
767 | Double_t bs=(binsizer+binsizez)/2.; |
f5f6da22 |
768 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
769 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); |
770 | if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i); |
70c95f95 |
771 | } |
f5f6da22 |
772 | fLines.Compress(); |
773 | AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast())); |
70c95f95 |
774 | |
f5f6da22 |
775 | if(fLines.GetEntriesFast()>1){ |
2ecf7a6f |
776 | retcode=0; |
27167524 |
777 | // find a first candidate for the primary vertex |
f5f6da22 |
778 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); |
70c95f95 |
779 | // make a further selection on tracklets based on this first candidate |
27167524 |
780 | fVert3D.GetXYZ(peak); |
781 | AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2])); |
f5f6da22 |
782 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
783 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); |
784 | if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i); |
70c95f95 |
785 | } |
f5f6da22 |
786 | fLines.Compress(); |
787 | AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast())); |
2ecf7a6f |
788 | if(fLines.GetEntriesFast()>1){// this new tracklet selection is used |
789 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); |
790 | } |
70c95f95 |
791 | } |
792 | return retcode; |
793 | } |
794 | |
6a0d56b8 |
795 | //________________________________________________________ |
0496a544 |
796 | void AliITSVertexer3D::SetMeanPPtSelTracks(Double_t fieldTesla){ |
5951029f |
797 | // Sets mean values of Pt based on the field |
798 | // for P (used in multiple scattering) the most probable value is used |
6a0d56b8 |
799 | if(TMath::Abs(fieldTesla-0.5)<0.01){ |
5951029f |
800 | SetMeanPSelTracks(0.375); |
6a0d56b8 |
801 | SetMeanPtSelTracks(0.630); |
802 | }else if(TMath::Abs(fieldTesla-0.4)<0.01){ |
5951029f |
803 | SetMeanPSelTracks(0.375); |
6a0d56b8 |
804 | SetMeanPtSelTracks(0.580); |
805 | }else if(TMath::Abs(fieldTesla-0.2)<0.01){ |
5951029f |
806 | SetMeanPSelTracks(0.375); |
6a0d56b8 |
807 | SetMeanPtSelTracks(0.530); |
808 | }else if(fieldTesla<0.00001){ |
5951029f |
809 | SetMeanPSelTracks(0.375); |
810 | SetMeanPtSelTracks(0.230); |
6a0d56b8 |
811 | }else{ |
812 | SetMeanPSelTracks(); |
813 | SetMeanPtSelTracks(); |
70c95f95 |
814 | } |
70c95f95 |
815 | } |
816 | |
7340864d |
817 | //________________________________________________________ |
818 | void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){ |
819 | // Finds bin with max contents in 3D histo of tracket intersections |
820 | TAxis *xax = histo->GetXaxis(); |
821 | TAxis *yax = histo->GetYaxis(); |
822 | TAxis *zax = histo->GetZaxis(); |
823 | peak[0]=0.; |
824 | peak[1]=0.; |
825 | peak[2]=0.; |
826 | nOfTracklets = 0; |
827 | nOfTimes=0; |
7d766abb |
828 | Int_t peakbin[3]={0,0,0}; |
829 | Int_t peak2bin[3]={-1,-1,-1}; |
830 | Int_t bc2=-1; |
7340864d |
831 | for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){ |
242c8025 |
832 | Double_t xval = xax->GetBinCenter(i); |
7340864d |
833 | for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){ |
242c8025 |
834 | Double_t yval = yax->GetBinCenter(j); |
7340864d |
835 | for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){ |
242c8025 |
836 | Double_t zval = zax->GetBinCenter(k); |
7340864d |
837 | Int_t bc =(Int_t)histo->GetBinContent(i,j,k); |
7d766abb |
838 | if(bc==0) continue; |
7340864d |
839 | if(bc>nOfTracklets){ |
7d766abb |
840 | nOfTracklets=bc; |
7340864d |
841 | peak[2] = zval; |
842 | peak[1] = yval; |
843 | peak[0] = xval; |
7d766abb |
844 | peakbin[2] = k; |
845 | peakbin[1] = j; |
846 | peakbin[0] = i; |
847 | peak2bin[2] = -1; |
848 | peak2bin[1] = -1; |
849 | peak2bin[0] = -1; |
850 | bc2=-1; |
7340864d |
851 | nOfTimes = 1; |
06a7cbee |
852 | }else if(bc==nOfTracklets){ |
7d766abb |
853 | if(TMath::Abs(i-peakbin[0])<=1 && TMath::Abs(j-peakbin[1])<=1 && TMath::Abs(k-peakbin[2])<=1){ |
854 | peak2bin[2] = k; |
855 | peak2bin[1] = j; |
856 | peak2bin[0] = i; |
857 | bc2=bc; |
858 | nOfTimes = 1; |
859 | }else{ |
860 | nOfTimes++; |
861 | } |
7340864d |
862 | } |
863 | } |
864 | } |
865 | } |
7d766abb |
866 | if(peak2bin[0]>=-1 && bc2!=-1){ // two contiguous peak-cells with same contents |
867 | peak[0]=0.5*(xax->GetBinCenter(peakbin[0])+xax->GetBinCenter(peak2bin[0])); |
868 | peak[1]=0.5*(yax->GetBinCenter(peakbin[1])+yax->GetBinCenter(peak2bin[1])); |
869 | peak[2]=0.5*(zax->GetBinCenter(peakbin[2])+zax->GetBinCenter(peak2bin[2])); |
870 | nOfTracklets+=bc2; |
871 | nOfTimes=1; |
872 | } |
7340864d |
873 | } |
cd1c2af1 |
874 | //________________________________________________________ |
875 | void AliITSVertexer3D::MarkUsedClusters(){ |
876 | // Mark clusters of tracklets used in vertex claulation |
877 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
878 | AliStrLine *lin = (AliStrLine*)fLines.At(i); |
879 | Int_t idClu1=lin->GetIdPoint(0); |
880 | Int_t idClu2=lin->GetIdPoint(1); |
881 | fUsedCluster.SetBitNumber(idClu1); |
882 | fUsedCluster.SetBitNumber(idClu2); |
883 | } |
884 | } |
885 | //________________________________________________________ |
886 | Int_t AliITSVertexer3D::RemoveTracklets(){ |
887 | // Remove trackelts close to first found vertex |
888 | Double_t vert[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()}; |
889 | Int_t nRemoved=0; |
890 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
891 | AliStrLine *lin = (AliStrLine*)fLines.At(i); |
892 | if(lin->GetDistFromPoint(vert)<fDCAforPileup){ |
893 | Int_t idClu1=lin->GetIdPoint(0); |
894 | Int_t idClu2=lin->GetIdPoint(1); |
895 | fUsedCluster.SetBitNumber(idClu1); |
896 | fUsedCluster.SetBitNumber(idClu2); |
897 | fLines.RemoveAt(i); |
898 | ++nRemoved; |
899 | } |
900 | } |
901 | fLines.Compress(); |
902 | return nRemoved; |
903 | } |
904 | //________________________________________________________ |
905 | void AliITSVertexer3D::FindOther3DVertices(TTree *itsClusterTree){ |
906 | // pileup identification based on 3D vertexing with not used clusters |
907 | MarkUsedClusters(); |
908 | fLines.Clear("C"); |
909 | Int_t nolines = FindTracklets(itsClusterTree,2); |
910 | if(nolines>=2){ |
911 | Int_t nr=RemoveTracklets(); |
912 | nolines-=nr; |
913 | if(nolines>=2){ |
914 | Int_t rc=Prepare3DVertex(2); |
915 | if(rc==0){ |
916 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); |
917 | if(fVert3D.GetNContributors()>=fMinTrackletsForPilup){ |
918 | fIsPileup=kTRUE; |
6b4d9537 |
919 | fNoVertices=2; |
920 | fVertArray = new AliESDVertex[2]; |
921 | fVertArray[0]=(*fCurrentVertex); |
922 | fVertArray[1]=fVert3D; |
cd1c2af1 |
923 | fZpuv=fVert3D.GetZv(); |
924 | fNTrpuv=fVert3D.GetNContributors(); |
925 | } |
926 | } |
927 | } |
928 | } |
929 | } |
930 | //______________________________________________________________________ |
931 | void AliITSVertexer3D::PileupFromZ(){ |
932 | // Calls the pileup algorithm of ALiITSVertexerZ |
933 | Int_t binmin, binmax; |
934 | Int_t nPeaks=AliITSVertexerZ::GetPeakRegion(fZHisto,binmin,binmax); |
935 | if(nPeaks==2)AliWarning("2 peaks found"); |
936 | Int_t firstPeakCont=0; |
242c8025 |
937 | Double_t firstPeakPos=0.; |
cd1c2af1 |
938 | for(Int_t i=binmin-1;i<=binmax+1;i++){ |
6b4d9537 |
939 | firstPeakCont+=static_cast<Int_t>(fZHisto->GetBinContent(i)); |
cd1c2af1 |
940 | firstPeakPos+=fZHisto->GetBinContent(i)*fZHisto->GetBinCenter(i); |
941 | } |
73315c4a |
942 | if(firstPeakCont>0){ |
943 | firstPeakPos/=firstPeakCont; |
944 | Int_t ncontr2=0; |
945 | if(firstPeakCont>fMinTrackletsForPilup){ |
946 | Float_t secPeakPos; |
947 | ncontr2=AliITSVertexerZ::FindSecondPeak(fZHisto,binmin,binmax,secPeakPos); |
948 | if(ncontr2>=fMinTrackletsForPilup){ |
949 | fIsPileup=kTRUE; |
6b4d9537 |
950 | fNoVertices=2; |
951 | AliESDVertex secondVert(secPeakPos,0.1,ncontr2); |
952 | fVertArray = new AliESDVertex[2]; |
953 | fVertArray[0]=(*fCurrentVertex); |
954 | fVertArray[1]=secondVert; |
73315c4a |
955 | fZpuv=secPeakPos; |
956 | fNTrpuv=ncontr2; |
957 | } |
cd1c2af1 |
958 | } |
959 | } |
960 | } |
70c95f95 |
961 | //________________________________________________________ |
962 | void AliITSVertexer3D::PrintStatus() const { |
963 | // Print current status |
87104b06 |
964 | printf("========= First step selections =====================\n"); |
6a0d56b8 |
965 | printf("Cut on diamond (Z) %f\n",fZCutDiamond); |
87104b06 |
966 | printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut); |
967 | printf("Loose cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut); |
6a0d56b8 |
968 | printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut); |
87104b06 |
969 | printf("========= Second step selections ====================\n"); |
970 | printf("Cut on tracklet-to-first-vertex Z distance %f\n",fMaxZCut); |
cd1c2af1 |
971 | printf("Max Phi difference: %f\n",fDiffPhiMax); |
87104b06 |
972 | printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut); |
973 | printf("Cut on tracklet DCA to beam axis (algo2) %f\n",fMaxRCut2); |
974 | printf("========= Pileup selections =========================\n"); |
cd1c2af1 |
975 | printf("Pileup algo: %d\n",fPileupAlgo); |
87104b06 |
976 | printf("Min DCA to 1st vertex for pileup (algo 0 and 1): %f\n",fDCAforPileup); |
977 | printf("Cut on distance between pair-vertices (algo 2): %f\n",fCutOnPairs); |
9801e044 |
978 | printf("Maximum number of clusters allowed on L1 or L2: %d\n",fMaxNumOfCl); |
6a0d56b8 |
979 | printf("=======================================================\n"); |
70c95f95 |
980 | } |