]>
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 | ||
efe6147d | 36 | const Int_t AliITSVertexer3D::fgkMaxNumOfClDefault = 300; |
37 | const Int_t AliITSVertexer3D::fgkMaxNumOfClRebinDefault = 500; | |
38 | const Int_t AliITSVertexer3D::fgkMaxNumOfClDownscaleDefault = 1000; | |
39 | const Float_t AliITSVertexer3D::fgk3DBinSizeDefault = 0.1; | |
9801e044 | 40 | |
70c95f95 | 41 | ClassImp(AliITSVertexer3D) |
42 | ||
6a0d56b8 | 43 | /* $Id$ */ |
44 | ||
70c95f95 | 45 | //______________________________________________________________________ |
4f40f207 | 46 | AliITSVertexer3D::AliITSVertexer3D(Double_t zcut): |
7d766abb | 47 | AliITSVertexer(), |
48 | fLines("AliStrLine",1000), | |
49 | fVert3D(), | |
50 | fCoarseDiffPhiCut(0.), | |
51 | fFineDiffPhiCut(0.), | |
52 | fCutOnPairs(0.), | |
53 | fCoarseMaxRCut(0.), | |
54 | fMaxRCut(0.), | |
87104b06 | 55 | fMaxRCut2(0.), |
7d766abb | 56 | fZCutDiamond(0.), |
57 | fMaxZCut(0.), | |
58 | fDCAcut(0.), | |
59 | fDiffPhiMax(0.), | |
60 | fMeanPSelTrk(0.), | |
61 | fMeanPtSelTrk(0.), | |
62 | fUsedCluster(kMaxCluPerMod*kNSPDMod), | |
63 | fZHisto(0), | |
64 | fDCAforPileup(0.), | |
c238ba33 | 65 | fDiffPhiforPileup(0.), |
7d766abb | 66 | fBinSizeR(0.), |
67 | fBinSizeZ(0.), | |
9801e044 | 68 | fPileupAlgo(0), |
507265fd | 69 | fMaxNumOfCl(fgkMaxNumOfClDefault), |
efe6147d | 70 | fMaxNumOfClForRebin(fgkMaxNumOfClRebinDefault), |
71 | fMaxNumOfClForDownScale(fgkMaxNumOfClDownscaleDefault), | |
72 | fNRecPLay1(0), | |
73 | fNRecPLay2(0), | |
74 | f3DBinSize(fgk3DBinSizeDefault), | |
507265fd | 75 | fDoDownScale(kFALSE), |
efe6147d | 76 | fGenerForDownScale(0), |
77 | f3DPeak(), | |
78 | fHighMultAlgo(1), | |
1d60ed44 | 79 | fSwitchAlgorithm(kFALSE), |
80 | fFallBack(kFALSE), | |
81 | fFallBackThreshold(0), | |
82 | fH3d(NULL), | |
83 | fH3dcs(NULL), | |
84 | fH3dfs(NULL), | |
85 | fH3dv(NULL) | |
6a0d56b8 | 86 | { |
70c95f95 | 87 | // Default constructor |
88 | SetCoarseDiffPhiCut(); | |
6b4d9537 | 89 | SetFineDiffPhiCut(); |
90 | SetCutOnPairs(); | |
05d1294c | 91 | SetCoarseMaxRCut(); |
70c95f95 | 92 | SetMaxRCut(); |
87104b06 | 93 | SetMaxRCutAlgo2(); |
4f40f207 | 94 | if(zcut>0.){ |
95 | SetZCutDiamond(zcut); | |
96 | } | |
97 | else { | |
98 | SetZCutDiamond(); | |
99 | } | |
05d1294c | 100 | SetMaxZCut(); |
7203e11a | 101 | SetDCACut(); |
70c95f95 | 102 | SetDiffPhiMax(); |
6a0d56b8 | 103 | SetMeanPSelTracks(); |
104 | SetMeanPtSelTracks(); | |
cd1c2af1 | 105 | SetMinDCAforPileup(); |
c238ba33 | 106 | SetDeltaPhiforPileup(); |
cd1c2af1 | 107 | SetPileupAlgo(); |
7d766abb | 108 | SetBinSizeR(); |
109 | SetBinSizeZ(); | |
507265fd | 110 | fGenerForDownScale=new TRandom3(987654321); |
70c95f95 | 111 | } |
112 | ||
f510fd70 | 113 | //______________________________________________________________________ |
114 | AliITSVertexer3D::AliITSVertexer3D(TRootIOCtor*): | |
115 | AliITSVertexer(), | |
116 | fLines("AliStrLine",1000), | |
117 | fVert3D(), | |
118 | fCoarseDiffPhiCut(0.), | |
119 | fFineDiffPhiCut(0.), | |
120 | fCutOnPairs(0.), | |
121 | fCoarseMaxRCut(0.), | |
122 | fMaxRCut(0.), | |
123 | fMaxRCut2(0.), | |
124 | fZCutDiamond(0.), | |
125 | fMaxZCut(0.), | |
126 | fDCAcut(0.), | |
127 | fDiffPhiMax(0.), | |
128 | fMeanPSelTrk(0.), | |
129 | fMeanPtSelTrk(0.), | |
130 | fUsedCluster(kMaxCluPerMod*kNSPDMod), | |
131 | fZHisto(0), | |
132 | fDCAforPileup(0.), | |
133 | fDiffPhiforPileup(0.), | |
134 | fBinSizeR(0.), | |
135 | fBinSizeZ(0.), | |
136 | fPileupAlgo(0), | |
137 | fMaxNumOfCl(fgkMaxNumOfClDefault), | |
138 | fMaxNumOfClForRebin(fgkMaxNumOfClRebinDefault), | |
139 | fMaxNumOfClForDownScale(fgkMaxNumOfClDownscaleDefault), | |
140 | fNRecPLay1(0), | |
141 | fNRecPLay2(0), | |
142 | f3DBinSize(fgk3DBinSizeDefault), | |
143 | fDoDownScale(kFALSE), | |
144 | fGenerForDownScale(0), | |
145 | f3DPeak(), | |
146 | fHighMultAlgo(1), | |
1d60ed44 | 147 | fSwitchAlgorithm(kFALSE), |
148 | fFallBack(kFALSE), | |
149 | fFallBackThreshold(0), | |
150 | fH3d(NULL), | |
151 | fH3dcs(NULL), | |
152 | fH3dfs(NULL), | |
153 | fH3dv(NULL) | |
f510fd70 | 154 | { |
155 | // I/O constructor | |
156 | ||
157 | ||
158 | } | |
159 | ||
70c95f95 | 160 | //______________________________________________________________________ |
161 | AliITSVertexer3D::~AliITSVertexer3D() { | |
162 | // Destructor | |
1d60ed44 | 163 | if(fH3d) delete fH3d; |
164 | if(fH3dcs) delete fH3dcs; | |
165 | if(fH3dfs) delete fH3dfs; | |
166 | if(fH3dv) delete fH3dv; | |
f5f6da22 | 167 | fLines.Clear("C"); |
cd1c2af1 | 168 | if(fZHisto) delete fZHisto; |
507265fd | 169 | if(fGenerForDownScale) delete fGenerForDownScale; |
70c95f95 | 170 | } |
171 | ||
27167524 | 172 | //______________________________________________________________________ |
173 | void AliITSVertexer3D::ResetVert3D(){ | |
f510fd70 | 174 | // Reset the fVert3D object and reset the used clusters |
6b4d9537 | 175 | ResetVertex(); |
27167524 | 176 | fVert3D.SetXv(0.); |
177 | fVert3D.SetYv(0.); | |
178 | fVert3D.SetZv(0.); | |
179 | fVert3D.SetDispersion(0.); | |
180 | fVert3D.SetNContributors(0); | |
cd1c2af1 | 181 | fUsedCluster.ResetAllBits(0); |
27167524 | 182 | } |
70c95f95 | 183 | //______________________________________________________________________ |
308c2f7c | 184 | AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree){ |
70c95f95 | 185 | // Defines the AliESDVertex for the current event |
1d60ed44 | 186 | |
187 | //cleanup | |
188 | if(fZHisto)fZHisto->Reset(); | |
27167524 | 189 | ResetVert3D(); |
308c2f7c | 190 | AliDebug(1,"FindVertexForCurrentEvent - 3D - PROCESSING NEXT EVENT"); |
6b4d9537 | 191 | fLines.Clear("C"); |
192 | fCurrentVertex = NULL; | |
70c95f95 | 193 | |
1d60ed44 | 194 | // fall back to VertexerZ if too many clusters on SPD first layer |
195 | AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance(); | |
196 | rpcont->FetchClusters(0,itsClusterTree); | |
197 | if(!rpcont->IsSPDActive()){ | |
198 | AliWarning("No SPD rec points found, 3D vertex not calculated"); | |
199 | return NULL; | |
200 | } | |
201 | Bool_t fallBack = kFALSE; // 3D algo , no fallback to vertexer Z | |
202 | if(fFallBack && (rpcont->GetNClustersInLayerFast(1)) > fFallBackThreshold) | |
203 | fallBack = kTRUE; | |
204 | if(!fallBack){ | |
205 | Int_t nolines = FindTracklets(itsClusterTree,0); | |
206 | Int_t rc; | |
207 | if(nolines>=2){ | |
208 | if(fSwitchAlgorithm) { | |
209 | rc = Prepare3DVertexPbPb(); | |
210 | FindVertex3D(itsClusterTree); | |
211 | } else { | |
212 | rc=Prepare3DVertex(0); | |
213 | if(fVert3D.GetNContributors()>0){ | |
214 | fLines.Clear("C"); | |
215 | nolines = FindTracklets(itsClusterTree,1); | |
216 | if(nolines>=2){ | |
217 | rc=Prepare3DVertex(1); | |
218 | if(fPileupAlgo == 2 && rc == 0) FindVertex3DIterative(); | |
219 | else if(fPileupAlgo!=2 && rc == 0) FindVertex3D(itsClusterTree); | |
220 | if(rc!=0) fVert3D.SetNContributors(0); // exclude this vertex | |
221 | } | |
efe6147d | 222 | } |
87104b06 | 223 | } |
224 | } | |
1d60ed44 | 225 | } // if(!fallBack) |
226 | if(fallBack || (!fCurrentVertex)){ | |
0b7bac55 | 227 | AliITSVertexerZ vertz(GetNominalPos()[0],GetNominalPos()[1]); |
b96ee725 | 228 | vertz.SetDetTypeRec(GetDetTypeRec()); |
0b7bac55 | 229 | AliDebug(1,"Call Vertexer Z\n"); |
cc282e54 | 230 | vertz.SetLowLimit(-fZCutDiamond); |
231 | vertz.SetHighLimit(fZCutDiamond); | |
0b7bac55 | 232 | AliESDVertex* vtxz = vertz.FindVertexForCurrentEvent(itsClusterTree); |
233 | if(vtxz){ | |
e690d4d0 | 234 | Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],vtxz->GetZ()}; |
0b7bac55 | 235 | Double_t covmatrix[6]; |
236 | vtxz->GetCovMatrix(covmatrix); | |
237 | Double_t chi2=99999.; | |
238 | Int_t nContr=vtxz->GetNContributors(); | |
239 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
bf2e0ad4 | 240 | fCurrentVertex->SetDispersion(vtxz->GetDispersion()); |
0b7bac55 | 241 | fCurrentVertex->SetTitle("vertexer: Z"); |
242 | fCurrentVertex->SetName("SPDVertexZ"); | |
243 | delete vtxz; | |
0b7bac55 | 244 | } |
245 | ||
1ff24d0a | 246 | } |
247 | if(fComputeMultiplicity) FindMultiplicity(itsClusterTree); | |
70c95f95 | 248 | return fCurrentVertex; |
249 | } | |
250 | ||
6b4d9537 | 251 | //______________________________________________________________________ |
252 | void AliITSVertexer3D::FindVertex3D(TTree *itsClusterTree){ | |
f510fd70 | 253 | // Instantiates the fCurrentVertex object. calle by FindVertexForCurrenEvent |
e690d4d0 | 254 | Double_t vRadius=TMath::Sqrt(fVert3D.GetX()*fVert3D.GetX()+fVert3D.GetY()*fVert3D.GetY()); |
6b4d9537 | 255 | if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){ |
e690d4d0 | 256 | Double_t position[3]={fVert3D.GetX(),fVert3D.GetY(),fVert3D.GetZ()}; |
6b4d9537 | 257 | Double_t covmatrix[6]; |
258 | fVert3D.GetCovMatrix(covmatrix); | |
259 | Double_t chi2=99999.; | |
260 | Int_t nContr=fVert3D.GetNContributors(); | |
261 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
262 | fCurrentVertex->SetTitle("vertexer: 3D"); | |
263 | fCurrentVertex->SetName("SPDVertex3D"); | |
264 | fCurrentVertex->SetDispersion(fVert3D.GetDispersion()); | |
265 | fNoVertices=1; | |
266 | ||
267 | switch(fPileupAlgo){ | |
268 | case 0: PileupFromZ(); break; | |
269 | case 1: FindOther3DVertices(itsClusterTree); break; | |
9ce18bbb | 270 | case 3: break; // no pileup algo |
6b4d9537 | 271 | default: AliError("Wrong pileup algorithm"); break; |
272 | } | |
273 | if(fNoVertices==1){ | |
049c07b3 | 274 | delete[] fVertArray; |
6b4d9537 | 275 | fVertArray = new AliESDVertex[1]; |
276 | fVertArray[0]=(*fCurrentVertex); | |
277 | } | |
278 | } | |
279 | } | |
280 | ||
281 | //______________________________________________________________________ | |
282 | void AliITSVertexer3D::FindVertex3DIterative(){ | |
f510fd70 | 283 | // find vertex if fPileupAlgo == 2 |
87104b06 | 284 | |
285 | Int_t nLines=fLines.GetEntriesFast(); | |
286 | Int_t maxPoints=nLines*(nLines-1)/2; | |
287 | Double_t* xP=new Double_t[maxPoints]; | |
288 | Double_t* yP=new Double_t[maxPoints]; | |
289 | Double_t* zP=new Double_t[maxPoints]; | |
290 | Int_t* index1=new Int_t[maxPoints]; | |
291 | Int_t* index2=new Int_t[maxPoints]; | |
e690d4d0 | 292 | Double_t xbeam=fVert3D.GetX(); |
293 | Double_t ybeam=fVert3D.GetY(); | |
87104b06 | 294 | |
295 | Int_t iPoint=0; | |
296 | for(Int_t ilin1=0; ilin1<nLines; ilin1++){ | |
297 | AliStrLine *l1 = (AliStrLine*)fLines.At(ilin1); | |
298 | for(Int_t ilin2=ilin1+1; ilin2<nLines; ilin2++){ | |
299 | AliStrLine *l2 = (AliStrLine*)fLines.At(ilin2); | |
300 | Double_t dca=l1->GetDCA(l2); | |
301 | if(dca > fDCAcut || dca<0.00001) continue; | |
302 | Double_t point[3]; | |
303 | Int_t retc = l1->Cross(l2,point); | |
304 | if(retc<0)continue; | |
305 | Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]); | |
306 | if(rad>fCoarseMaxRCut)continue; | |
307 | Double_t distFromBeam=TMath::Sqrt((point[0]-xbeam)*(point[0]-xbeam)+(point[1]-ybeam)*(point[1]-ybeam)); | |
308 | if(distFromBeam>fMaxRCut2) continue; | |
309 | xP[iPoint]=point[0]; | |
310 | yP[iPoint]=point[1]; | |
311 | zP[iPoint]=point[2]; | |
312 | index1[iPoint]=ilin1; | |
313 | index2[iPoint]=ilin2; | |
314 | iPoint++; | |
315 | } | |
316 | } | |
317 | Int_t npoints=iPoint++; | |
318 | Int_t index=0; | |
319 | Short_t* mask=new Short_t[npoints]; | |
320 | for(Int_t ip=0;ip<npoints;ip++) mask[ip]=-1; | |
321 | ||
322 | for(Int_t ip1=0;ip1<npoints;ip1++){ | |
323 | if(mask[ip1]==-1) mask[ip1]=index++; | |
324 | for(Int_t ip2=ip1+1; ip2<npoints; ip2++){ | |
325 | if(mask[ip2]==mask[ip1] && mask[ip2]!=-1) continue; | |
326 | Double_t dist2=(xP[ip1]-xP[ip2])*(xP[ip1]-xP[ip2]); | |
327 | dist2+=(yP[ip1]-yP[ip2])*(yP[ip1]-yP[ip2]); | |
328 | dist2+=(zP[ip1]-zP[ip2])*(zP[ip1]-zP[ip2]); | |
329 | if(dist2<fCutOnPairs*fCutOnPairs){ | |
330 | if(mask[ip2]==-1) mask[ip2]=mask[ip1]; | |
331 | else{ | |
332 | for(Int_t ip=0; ip<npoints;ip++){ | |
333 | if(mask[ip]==mask[ip2]) mask[ip]=mask[ip1]; | |
334 | } | |
335 | } | |
336 | } | |
337 | } | |
338 | } | |
339 | ||
340 | ||
341 | // Count multiplicity of trackelts in clusters | |
342 | UInt_t* isIndUsed=new UInt_t[index+1]; | |
343 | for(Int_t ind=0;ind<index+1;ind++) isIndUsed[ind]=0; | |
344 | for(Int_t ip=0; ip<npoints;ip++){ | |
345 | Int_t ind=mask[ip]; | |
346 | isIndUsed[ind]++; | |
347 | } | |
348 | ||
349 | // Count clusters/vertices and sort according to multiplicity | |
350 | Int_t nClusters=0; | |
351 | Int_t* sortedIndex=new Int_t[index+1]; | |
352 | for(Int_t ind1=0;ind1<index+1;ind1++){ | |
353 | if(isIndUsed[ind1]<=1) isIndUsed[ind1]=0; | |
354 | else nClusters++; | |
355 | UInt_t cap=9999999; | |
356 | if(ind1>0) cap=isIndUsed[sortedIndex[ind1-1]]; | |
357 | UInt_t bigger=0; | |
358 | Int_t biggerindex=-1; | |
359 | for(Int_t ind2=0;ind2<index+1;ind2++){ | |
360 | Bool_t use=kTRUE; | |
361 | for(Int_t ind3=0; ind3<ind1; ind3++) | |
362 | if(ind2==sortedIndex[ind3]) use=kFALSE; | |
363 | if(use && isIndUsed[ind2]>bigger && isIndUsed[ind2]<=cap){ | |
364 | bigger=isIndUsed[ind2]; | |
365 | biggerindex=ind2; | |
366 | } | |
367 | } | |
368 | sortedIndex[ind1]=biggerindex; | |
369 | } | |
370 | AliDebug(3,Form("Number of clusters before merging = %d\n",nClusters)); | |
371 | ||
372 | // Assign lines to clusters/vertices and merge clusters which share 1 line | |
373 | Int_t nClustersAfterMerge=nClusters; | |
374 | Int_t* belongsTo=new Int_t[nLines]; | |
375 | for(Int_t ilin=0; ilin<nLines; ilin++) belongsTo[ilin]=-1; | |
376 | for(Int_t iclu=0;iclu<nClusters;iclu++){ | |
377 | Int_t actualCluIndex=iclu; | |
378 | for(Int_t ip=0; ip<npoints;ip++){ | |
379 | if(mask[ip]==sortedIndex[iclu]){ | |
380 | Int_t ind1=index1[ip]; | |
381 | if(belongsTo[ind1]==-1) belongsTo[ind1]=actualCluIndex; | |
382 | else if(belongsTo[ind1]<actualCluIndex){ | |
383 | Int_t newCluIndex=belongsTo[ind1]; | |
384 | for(Int_t ilin=0; ilin<nLines; ilin++){ | |
385 | if(belongsTo[ilin]==actualCluIndex) belongsTo[ilin]=newCluIndex; | |
386 | } | |
387 | AliDebug(10,Form("Merged cluster %d with %d\n",actualCluIndex,newCluIndex)); | |
388 | actualCluIndex=newCluIndex; | |
389 | nClustersAfterMerge--; | |
390 | } | |
391 | Int_t ind2=index2[ip]; | |
392 | if(belongsTo[ind2]==-1) belongsTo[ind2]=actualCluIndex; | |
393 | else if(belongsTo[ind2]<actualCluIndex){ | |
394 | Int_t newCluIndex=belongsTo[ind2]; | |
395 | for(Int_t ilin=0; ilin<nLines; ilin++){ | |
396 | if(belongsTo[ilin]==actualCluIndex) belongsTo[ilin]=newCluIndex; | |
397 | } | |
398 | AliDebug(10,Form("Merged cluster %d with %d\n",actualCluIndex,newCluIndex)); | |
399 | actualCluIndex=newCluIndex; | |
400 | nClustersAfterMerge--; | |
401 | } | |
402 | } | |
403 | } | |
404 | } | |
405 | AliDebug(3,Form("Number of clusters after merging = %d\n",nClustersAfterMerge)); | |
406 | ||
407 | // Count lines associated to each cluster/vertex | |
408 | UInt_t *cluSize=new UInt_t[nClusters]; | |
409 | for(Int_t iclu=0;iclu<nClusters;iclu++){ | |
410 | cluSize[iclu]=0; | |
411 | for(Int_t ilin=0; ilin<nLines; ilin++){ | |
412 | if(belongsTo[ilin]==iclu) cluSize[iclu]++; | |
413 | } | |
414 | } | |
415 | ||
416 | // Count good vertices (>1 associated tracklet) | |
417 | UInt_t nGoodVert=0; | |
418 | for(Int_t iclu=0;iclu<nClusters;iclu++){ | |
419 | AliDebug(3,Form("Vertex %d Size=%d\n",iclu,cluSize[iclu])); | |
420 | if(cluSize[iclu]>1) nGoodVert++; | |
421 | } | |
422 | ||
423 | AliDebug(1,Form("Number of good vertices = %d\n",nGoodVert)); | |
424 | // Calculate vertex coordinates for each cluster | |
425 | if(nGoodVert>0){ | |
426 | fVertArray = new AliESDVertex[nGoodVert]; | |
427 | Int_t iVert=0; | |
428 | for(Int_t iclu=0;iclu<nClusters;iclu++){ | |
429 | Int_t size=cluSize[iclu]; | |
430 | if(size>1){ | |
431 | AliStrLine **arrlin = new AliStrLine*[size]; | |
432 | Int_t nFilled=0; | |
433 | for(Int_t ilin=0; ilin<nLines; ilin++){ | |
434 | if(belongsTo[ilin]==iclu){ | |
435 | arrlin[nFilled++] = dynamic_cast<AliStrLine*>(fLines[ilin]); | |
436 | } | |
437 | } | |
438 | AliDebug(3,Form("Vertex %d N associated tracklets = %d out of %d\n",iVert,size,nFilled)); | |
439 | ||
440 | fVertArray[iVert]=AliVertexerTracks::TrackletVertexFinder(arrlin,nFilled); | |
441 | Double_t peak[3]; | |
442 | fVertArray[iVert].GetXYZ(peak); | |
443 | AliStrLine **arrlin2 = new AliStrLine*[size]; | |
444 | Int_t nFilled2=0; | |
445 | for(Int_t i=0; i<nFilled;i++){ | |
446 | AliStrLine *l1 = arrlin[i]; | |
447 | if(l1->GetDistFromPoint(peak)< fDCAcut) | |
448 | arrlin2[nFilled2++] = dynamic_cast<AliStrLine*>(l1); | |
449 | } | |
450 | if(nFilled2>1){ | |
451 | AliDebug(3,Form("Vertex %d recalculated with %d tracklets\n",iVert,nFilled2)); | |
452 | fVertArray[iVert]=AliVertexerTracks::TrackletVertexFinder(arrlin2,nFilled2); | |
453 | } | |
454 | delete [] arrlin; | |
455 | delete [] arrlin2; | |
456 | ++iVert; | |
457 | } | |
458 | } | |
459 | ||
460 | if(nGoodVert > 1){ | |
461 | fIsPileup = kTRUE; | |
462 | fNTrpuv = fVertArray[1].GetNContributors(); | |
e690d4d0 | 463 | fZpuv = fVertArray[1].GetZ(); |
87104b06 | 464 | } |
465 | ||
e690d4d0 | 466 | Double_t vRadius=TMath::Sqrt(fVertArray[0].GetX()*fVertArray[0].GetX()+fVertArray[0].GetY()*fVertArray[0].GetY()); |
87104b06 | 467 | if(vRadius<GetPipeRadius() && fVertArray[0].GetNContributors()>0){ |
e690d4d0 | 468 | Double_t position[3]={fVertArray[0].GetX(),fVertArray[0].GetY(),fVertArray[0].GetZ()}; |
87104b06 | 469 | Double_t covmatrix[6]; |
470 | fVertArray[0].GetCovMatrix(covmatrix); | |
471 | Double_t chi2=99999.; | |
472 | Int_t nContr=fVertArray[0].GetNContributors(); | |
473 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
474 | fCurrentVertex->SetTitle("vertexer: 3D"); | |
475 | fCurrentVertex->SetName("SPDVertex3D"); | |
476 | fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion()); | |
477 | } | |
478 | } | |
479 | ||
480 | delete [] index1; | |
481 | delete [] index2; | |
482 | delete [] mask; | |
483 | delete [] isIndUsed; | |
484 | delete [] sortedIndex; | |
485 | delete [] belongsTo; | |
486 | delete [] cluSize; | |
487 | delete [] xP; | |
488 | delete [] yP; | |
489 | delete [] zP; | |
490 | } | |
491 | //______________________________________________________________________ | |
492 | void AliITSVertexer3D::FindVertex3DIterativeMM(){ | |
6b4d9537 | 493 | // Defines the AliESDVertex for the current event |
494 | Int_t numsor=fLines.GetEntriesFast()*(fLines.GetEntriesFast()-1)/2; | |
495 | //cout<<"AliITSVertexer3D::FindVertexForCurentEvent: Number of tracklets selected for vertexing "<<fLines.GetEntriesFast()<<"; Number of pairs: "<<numsor<<endl; | |
496 | AliITSSortTrkl srt(fLines,numsor,fCutOnPairs,fCoarseMaxRCut); | |
497 | srt.FindClusters(); | |
498 | AliInfo(Form("Number of vertices: %d",srt.GetNumberOfClusters())); | |
499 | ||
500 | fNoVertices = srt.GetNumberOfClusters(); | |
501 | //printf("fNoVertices = %d \n",fNoVertices); | |
502 | if(fNoVertices>0){ | |
503 | fVertArray = new AliESDVertex[fNoVertices]; | |
504 | for(Int_t kk=0; kk<srt.GetNumberOfClusters(); kk++){ | |
505 | Int_t size = 0; | |
506 | Int_t *labels = srt.GetTrackletsLab(kk,size); | |
507 | /* | |
508 | Int_t *pairs = srt.GetClusters(kk); | |
509 | Int_t nopai = srt.GetSizeOfCluster(kk); | |
510 | cout<<"***** Vertex number "<<kk<<". Pairs: \n"; | |
511 | for(Int_t jj=0;jj<nopai;jj++){ | |
512 | cout<<pairs[jj]<<" - "; | |
513 | if(jj>0 & jj%8==0)cout<<endl; | |
514 | } | |
515 | cout<<endl; | |
516 | cout<<"***** Vertex number "<<kk<<". Labels: \n"; | |
517 | */ | |
518 | AliStrLine **tclo = new AliStrLine* [size]; | |
519 | for(Int_t jj=0;jj<size;jj++){ | |
520 | // cout<<labels[jj]<<" - "; | |
521 | // if(jj>0 & jj%8==0)cout<<endl; | |
522 | tclo[jj] = dynamic_cast<AliStrLine*>(fLines[labels[jj]]); | |
523 | } | |
524 | // cout<<endl; | |
525 | delete []labels; | |
526 | fVertArray[kk]=AliVertexerTracks::TrackletVertexFinder(tclo,size); | |
527 | delete [] tclo; | |
528 | // fVertArray[kk].PrintStatus(); | |
529 | if(kk == 1){ | |
530 | // at least one second vertex is present | |
531 | fIsPileup = kTRUE; | |
532 | fNTrpuv = fVertArray[kk].GetNContributors(); | |
e690d4d0 | 533 | fZpuv = fVertArray[kk].GetZ(); |
6b4d9537 | 534 | } |
535 | } | |
e690d4d0 | 536 | Double_t vRadius=TMath::Sqrt(fVertArray[0].GetX()*fVertArray[0].GetX()+fVertArray[0].GetY()*fVertArray[0].GetY()); |
6b4d9537 | 537 | if(vRadius<GetPipeRadius() && fVertArray[0].GetNContributors()>0){ |
e690d4d0 | 538 | Double_t position[3]={fVertArray[0].GetX(),fVertArray[0].GetY(),fVertArray[0].GetZ()}; |
6b4d9537 | 539 | Double_t covmatrix[6]; |
540 | fVertArray[0].GetCovMatrix(covmatrix); | |
541 | Double_t chi2=99999.; | |
542 | Int_t nContr=fVertArray[0].GetNContributors(); | |
543 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
544 | fCurrentVertex->SetTitle("vertexer: 3D"); | |
545 | fCurrentVertex->SetName("SPDVertex3D"); | |
546 | fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion()); | |
547 | } | |
548 | } | |
549 | ||
550 | } | |
551 | ||
552 | //______________________________________________________________________ | |
553 | Bool_t AliITSVertexer3D::DistBetweenVertices(AliESDVertex &a, AliESDVertex &b, Double_t test, Double_t &dist){ | |
554 | // method to compare the distance between vertices a and b with "test" | |
555 | //it returns kTRUE is the distance is less or equal to test | |
556 | dist = (a.GetX()-b.GetX()) * (a.GetX()-b.GetX()); | |
557 | dist += (a.GetY()-b.GetY()) * (a.GetY()-b.GetY()); | |
558 | dist += (a.GetZ()-b.GetZ()) * (a.GetZ()-b.GetZ()); | |
559 | dist = TMath::Sqrt(dist); | |
560 | if(dist <= test)return kTRUE; | |
561 | return kFALSE; | |
562 | } | |
563 | ||
564 | ||
70c95f95 | 565 | //______________________________________________________________________ |
308c2f7c | 566 | Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){ |
70c95f95 | 567 | // All the possible combinations between recpoints on layer 1and 2 are |
568 | // considered. Straight lines (=tracklets)are formed. | |
569 | // The tracklets are processed in Prepare3DVertex | |
70c95f95 | 570 | |
70c95f95 | 571 | TClonesArray *itsRec = 0; |
cd1c2af1 | 572 | if(optCuts==0) fZHisto->Reset(); |
efe6147d | 573 | // gc1 are local and global coordinates for layer 1 |
0496a544 | 574 | Float_t gc1f[3]={0.,0.,0.}; |
575 | Double_t gc1[3]={0.,0.,0.}; | |
6b4d9537 | 576 | // gc2 are local and global coordinates for layer 2 |
0496a544 | 577 | Float_t gc2f[3]={0.,0.,0.}; |
578 | Double_t gc2[3]={0.,0.,0.}; | |
b21c1af0 | 579 | AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance(); |
3457dd2d | 580 | rpcont->FetchClusters(0,itsClusterTree); |
b21c1af0 | 581 | if(!rpcont->IsSPDActive()){ |
582 | AliWarning("No SPD rec points found, 3D vertex not calculated"); | |
6b4d9537 | 583 | return -1; |
584 | } | |
70c95f95 | 585 | |
586 | // Set values for cuts | |
0496a544 | 587 | Double_t xbeam=GetNominalPos()[0]; |
588 | Double_t ybeam=GetNominalPos()[1]; | |
589 | Double_t zvert=0.; | |
590 | Double_t deltaPhi=fCoarseDiffPhiCut; | |
591 | Double_t deltaR=fCoarseMaxRCut; | |
592 | Double_t dZmax=fZCutDiamond; | |
87104b06 | 593 | if(optCuts==1){ |
e690d4d0 | 594 | xbeam=fVert3D.GetX(); |
595 | ybeam=fVert3D.GetY(); | |
596 | zvert=fVert3D.GetZ(); | |
70c95f95 | 597 | deltaPhi = fDiffPhiMax; |
05d1294c | 598 | deltaR=fMaxRCut; |
599 | dZmax=fMaxZCut; | |
87104b06 | 600 | if(fPileupAlgo == 2){ |
601 | dZmax=fZCutDiamond; | |
602 | deltaR=fMaxRCut2; | |
603 | } | |
6b4d9537 | 604 | } else if(optCuts==2){ |
e690d4d0 | 605 | xbeam=fVert3D.GetX(); |
606 | ybeam=fVert3D.GetY(); | |
c238ba33 | 607 | deltaPhi = fDiffPhiforPileup; |
cd1c2af1 | 608 | deltaR=fMaxRCut; |
70c95f95 | 609 | } |
cd1c2af1 | 610 | |
efe6147d | 611 | fNRecPLay1=rpcont->GetNClustersInLayerFast(1); |
612 | fNRecPLay2=rpcont->GetNClustersInLayerFast(2); | |
613 | if(fNRecPLay1 == 0 || fNRecPLay2 == 0){ | |
614 | AliDebug(1,Form("No RecPoints in at least one SPD layer (%d %d)",fNRecPLay1,fNRecPLay2)); | |
70c95f95 | 615 | return -1; |
616 | } | |
efe6147d | 617 | AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",fNRecPLay1,fNRecPLay2)); |
507265fd | 618 | fDoDownScale=kFALSE; |
efe6147d | 619 | fSwitchAlgorithm=kFALSE; |
620 | ||
507265fd | 621 | Float_t factDownScal=1.; |
622 | Int_t origLaddersOnLayer2=fLadOnLay2; | |
623 | ||
efe6147d | 624 | switch(fHighMultAlgo){ |
625 | case 0: | |
626 | if(fNRecPLay1>fMaxNumOfClForDownScale || fNRecPLay2>fMaxNumOfClForDownScale){ | |
627 | if(optCuts==2) return -1; // do not try to search for pileup | |
628 | SetLaddersOnLayer2(2); | |
629 | fDoDownScale=kTRUE; | |
630 | factDownScal=(Float_t)fMaxNumOfClForDownScale*(Float_t)fMaxNumOfClForDownScale/(Float_t)fNRecPLay1/(Float_t)fNRecPLay2; | |
631 | if(optCuts==1){ | |
632 | factDownScal*=(fCoarseDiffPhiCut/fDiffPhiMax)*10; | |
633 | if(factDownScal>1.){ | |
634 | fDoDownScale=kFALSE; | |
635 | SetLaddersOnLayer2(origLaddersOnLayer2); | |
636 | } | |
507265fd | 637 | } |
efe6147d | 638 | if(fDoDownScale)AliDebug(1,Form("Too many recpoints on SPD(%d %d ), downscale by %f",fNRecPLay1,fNRecPLay2,factDownScal)); |
639 | } | |
640 | break; | |
641 | case 1: | |
642 | if(fNRecPLay1>fMaxNumOfCl || fNRecPLay2>fMaxNumOfCl) { | |
643 | if(optCuts==2) return -1; // do not try to search for pileup | |
644 | fSwitchAlgorithm=kTRUE; | |
645 | } | |
646 | break; | |
647 | default: break; // no pileup algo | |
648 | } | |
649 | ||
650 | if(!fDoDownScale && !fSwitchAlgorithm){ | |
651 | if(fNRecPLay1>fMaxNumOfClForRebin || fNRecPLay2>fMaxNumOfClForRebin){ | |
652 | SetLaddersOnLayer2(2); | |
507265fd | 653 | } |
9801e044 | 654 | } |
70c95f95 | 655 | |
656 | Double_t a[3]={xbeam,ybeam,0.}; | |
657 | Double_t b[3]={xbeam,ybeam,10.}; | |
658 | AliStrLine zeta(a,b,kTRUE); | |
322e007a | 659 | static Double_t bField=TMath::Abs(AliTracker::GetBz()/10.); //T |
6a0d56b8 | 660 | SetMeanPPtSelTracks(bField); |
87104b06 | 661 | |
70c95f95 | 662 | Int_t nolines = 0; |
663 | // Loop on modules of layer 1 | |
32f7748d | 664 | Int_t firstL1 = TMath::Max(0,AliITSgeomTGeo::GetModuleIndex(1,1,1)); |
1cc75a0b | 665 | Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1; |
6a0d56b8 | 666 | for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1 |
8c42830a | 667 | if(!fUseModule[modul1]) continue; |
507265fd | 668 | |
410cfb90 | 669 | UShort_t ladder=modul1/4+1; // ladders are numbered starting from 1 |
b21c1af0 | 670 | TClonesArray *prpl1=rpcont->UncheckedGetClusters(modul1); |
671 | Int_t nrecp1 = prpl1->GetEntries(); | |
70c95f95 | 672 | for(Int_t j=0;j<nrecp1;j++){ |
cd1c2af1 | 673 | if(j>kMaxCluPerMod) continue; |
674 | UShort_t idClu1=modul1*kMaxCluPerMod+j; | |
675 | if(fUsedCluster.TestBitNumber(idClu1)) continue; | |
efe6147d | 676 | if(fDoDownScale && !fSwitchAlgorithm){ |
677 | if(fGenerForDownScale->Rndm()>factDownScal) continue; | |
678 | } | |
b21c1af0 | 679 | AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1->At(j); |
0496a544 | 680 | recp1->GetGlobalXYZ(gc1f); |
681 | for(Int_t ico=0;ico<3;ico++)gc1[ico]=gc1f[ico]; | |
682 | ||
6a0d56b8 | 683 | Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam); |
70c95f95 | 684 | if(phi1<0)phi1=2*TMath::Pi()+phi1; |
27167524 | 685 | for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){ |
686 | for(Int_t k=0;k<4;k++){ | |
687 | Int_t ladmod=fLadders[ladder-1]+ladl2; | |
32e63e47 | 688 | if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2); |
689 | Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1); | |
e917bef4 | 690 | if(modul2<0)continue; |
8c42830a | 691 | if(!fUseModule[modul2]) continue; |
b21c1af0 | 692 | itsRec=rpcont->UncheckedGetClusters(modul2); |
27167524 | 693 | Int_t nrecp2 = itsRec->GetEntries(); |
694 | for(Int_t j2=0;j2<nrecp2;j2++){ | |
cd1c2af1 | 695 | if(j2>kMaxCluPerMod) continue; |
696 | UShort_t idClu2=modul2*kMaxCluPerMod+j2; | |
697 | if(fUsedCluster.TestBitNumber(idClu2)) continue; | |
507265fd | 698 | |
6a0d56b8 | 699 | AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2); |
0496a544 | 700 | recp2->GetGlobalXYZ(gc2f); |
701 | for(Int_t ico=0;ico<3;ico++)gc2[ico]=gc2f[ico]; | |
27167524 | 702 | Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam); |
703 | if(phi2<0)phi2=2*TMath::Pi()+phi2; | |
704 | Double_t diff = TMath::Abs(phi2-phi1); | |
705 | if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; | |
54091147 | 706 | if(optCuts==0 && diff<fDiffPhiforPileup){ |
0496a544 | 707 | Double_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); |
708 | Double_t zc1=gc1[2]; | |
709 | Double_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); | |
710 | Double_t zc2=gc2[2]; | |
711 | Double_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius | |
cd1c2af1 | 712 | fZHisto->Fill(zr0); |
713 | } | |
27167524 | 714 | if(diff>deltaPhi)continue; |
6a0d56b8 | 715 | AliStrLine line(gc1,gc2,kTRUE); |
27167524 | 716 | Double_t cp[3]; |
717 | Int_t retcode = line.Cross(&zeta,cp); | |
718 | if(retcode<0)continue; | |
719 | Double_t dca = line.GetDCA(&zeta); | |
720 | if(dca<0.) continue; | |
721 | if(dca>deltaR)continue; | |
722 | Double_t deltaZ=cp[2]-zvert; | |
723 | if(TMath::Abs(deltaZ)>dZmax)continue; | |
6a0d56b8 | 724 | |
507265fd | 725 | |
6a0d56b8 | 726 | if(nolines == 0){ |
f5f6da22 | 727 | if(fLines.GetEntriesFast()>0)fLines.Clear("C"); |
6a0d56b8 | 728 | } |
729 | Float_t cov[6]; | |
730 | recp2->GetGlobalCov(cov); | |
731 | ||
cd1c2af1 | 732 | |
0496a544 | 733 | Double_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); |
734 | Double_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); | |
735 | Double_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction | |
6a0d56b8 | 736 | |
0496a544 | 737 | Double_t curvErr=0; |
6a0d56b8 | 738 | if(bField>0.00001){ |
0496a544 | 739 | Double_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm |
440198c1 | 740 | Double_t dRad=TMath::Sqrt((gc1[0]-gc2[0])*(gc1[0]-gc2[0])+(gc1[1]-gc2[1])*(gc1[1]-gc2[1])); |
0496a544 | 741 | Double_t aux=dRad/2.+rad1; |
6a0d56b8 | 742 | curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm |
743 | } | |
0496a544 | 744 | Double_t sigmasq[3]; |
6a0d56b8 | 745 | sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor; |
746 | sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor; | |
747 | sigmasq[2]=cov[5]*factor*factor; | |
748 | ||
749 | // Multiple scattering | |
0496a544 | 750 | Double_t pOverMass=fMeanPSelTrk/0.140; |
751 | Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass); | |
752 | Double_t p2=fMeanPSelTrk*fMeanPSelTrk; | |
753 | Double_t rBP=GetPipeRadius(); | |
754 | Double_t dBP=0.08/35.3; // 800 um of Be | |
755 | Double_t dL1=0.01; //approx. 1% of radiation length | |
756 | Double_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP; | |
757 | Double_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1; | |
758 | Double_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1)); | |
759 | Double_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP)); | |
6a0d56b8 | 760 | for(Int_t ico=0; ico<3;ico++){ |
5951029f | 761 | sigmasq[ico]+=rtantheta1*rtantheta1*factor*factor/3.; |
762 | sigmasq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.; | |
6a0d56b8 | 763 | } |
0496a544 | 764 | Double_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.}; |
6a0d56b8 | 765 | if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0]; |
766 | if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1]; | |
767 | if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2]; | |
cd1c2af1 | 768 | new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE,idClu1,idClu2); |
6a0d56b8 | 769 | |
27167524 | 770 | } |
70c95f95 | 771 | } |
70c95f95 | 772 | } |
773 | } | |
70c95f95 | 774 | } |
507265fd | 775 | |
efe6147d | 776 | SetLaddersOnLayer2(origLaddersOnLayer2); |
507265fd | 777 | |
70c95f95 | 778 | if(nolines == 0)return -2; |
779 | return nolines; | |
780 | } | |
781 | ||
70c95f95 | 782 | //______________________________________________________________________ |
05d1294c | 783 | Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ |
70c95f95 | 784 | // Finds the 3D vertex information using tracklets |
785 | Int_t retcode = -1; | |
0496a544 | 786 | Double_t xbeam=GetNominalPos()[0]; |
787 | Double_t ybeam=GetNominalPos()[1]; | |
788 | Double_t zvert=0.; | |
789 | Double_t deltaR=fCoarseMaxRCut; | |
0496a544 | 790 | Double_t dZmax=fZCutDiamond; |
cd1c2af1 | 791 | if(optCuts==1){ |
e690d4d0 | 792 | xbeam=fVert3D.GetX(); |
793 | ybeam=fVert3D.GetY(); | |
794 | zvert=fVert3D.GetZ(); | |
05d1294c | 795 | deltaR=fMaxRCut; |
796 | dZmax=fMaxZCut; | |
87104b06 | 797 | if(fPileupAlgo == 2){ |
798 | dZmax=fZCutDiamond; | |
799 | deltaR=fMaxRCut2; | |
800 | } | |
cd1c2af1 | 801 | }else if(optCuts==2){ |
e690d4d0 | 802 | xbeam=fVert3D.GetX(); |
803 | ybeam=fVert3D.GetY(); | |
cd1c2af1 | 804 | deltaR=fMaxRCut; |
05d1294c | 805 | } |
806 | ||
507265fd | 807 | Double_t origBinSizeR=fBinSizeR; |
808 | Double_t origBinSizeZ=fBinSizeZ; | |
efe6147d | 809 | Bool_t rebinned=kFALSE; |
507265fd | 810 | if(fDoDownScale){ |
811 | SetBinSizeR(0.05); | |
812 | SetBinSizeZ(0.05); | |
efe6147d | 813 | rebinned=kTRUE; |
814 | }else{ | |
815 | if(optCuts==0 && (fNRecPLay1>fMaxNumOfClForRebin || fNRecPLay2>fMaxNumOfClForRebin)){ | |
816 | SetBinSizeR(0.1); | |
817 | SetBinSizeZ(0.2); | |
818 | rebinned=kTRUE; | |
819 | } | |
507265fd | 820 | } |
0496a544 | 821 | Double_t rl=-fCoarseMaxRCut; |
822 | Double_t rh=fCoarseMaxRCut; | |
823 | Double_t zl=-fZCutDiamond; | |
824 | Double_t zh=fZCutDiamond; | |
7d766abb | 825 | Int_t nbr=(Int_t)((rh-rl)/fBinSizeR+0.0001); |
826 | Int_t nbz=(Int_t)((zh-zl)/fBinSizeZ+0.0001); | |
827 | Int_t nbrcs=(Int_t)((rh-rl)/(fBinSizeR*2.)+0.0001); | |
828 | Int_t nbzcs=(Int_t)((zh-zl)/(fBinSizeZ*2.)+0.0001); | |
1d60ed44 | 829 | if(!fH3d){ |
830 | fH3d = new TH3F("fH3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh); | |
831 | fH3d->SetDirectory(0); | |
832 | }else{ | |
833 | fH3d->SetBins(nbr,rl,rh,nbr,rl,rh,nbz,zl,zh); | |
834 | } | |
835 | if(!fH3dcs){ | |
836 | fH3dcs = new TH3F("fH3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh); | |
837 | fH3dcs->SetDirectory(0); | |
838 | }else{ | |
839 | fH3dcs->SetBins(nbr,rl,rh,nbr,rl,rh,nbz,zl,zh); | |
840 | } | |
87104b06 | 841 | |
70c95f95 | 842 | // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed) |
3457dd2d | 843 | Int_t vsiz = fLines.GetEntriesFast(); |
844 | Int_t *validate = new Int_t [vsiz]; | |
845 | for(Int_t i=0; i<vsiz;i++)validate[i]=0; | |
846 | for(Int_t i=0; i<vsiz-1;i++){ | |
f5f6da22 | 847 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); |
848 | for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){ | |
849 | AliStrLine *l2 = (AliStrLine*)fLines.At(j); | |
0496a544 | 850 | Double_t dca=l1->GetDCA(l2); |
05d1294c | 851 | if(dca > fDCAcut || dca<0.00001) continue; |
70c95f95 | 852 | Double_t point[3]; |
853 | Int_t retc = l1->Cross(l2,point); | |
854 | if(retc<0)continue; | |
6a0d56b8 | 855 | Double_t deltaZ=point[2]-zvert; |
efe6147d | 856 | if(TMath::Abs(deltaZ)>dZmax)continue; |
70c95f95 | 857 | Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]); |
05d1294c | 858 | if(rad>fCoarseMaxRCut)continue; |
859 | Double_t deltaX=point[0]-xbeam; | |
860 | Double_t deltaY=point[1]-ybeam; | |
05d1294c | 861 | Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY); |
05d1294c | 862 | if(raddist>deltaR)continue; |
27167524 | 863 | validate[i]=1; |
864 | validate[j]=1; | |
1d60ed44 | 865 | fH3d->Fill(point[0],point[1],point[2]); |
866 | fH3dcs->Fill(point[0],point[1],point[2]); | |
70c95f95 | 867 | } |
868 | } | |
869 | ||
70c95f95 | 870 | Int_t numbtracklets=0; |
3457dd2d | 871 | for(Int_t i=0; i<vsiz;i++)if(validate[i]>=1)numbtracklets++; |
6b4d9537 | 872 | if(numbtracklets<2){ |
873 | delete [] validate; | |
1d60ed44 | 874 | fH3d->Reset(); |
875 | fH3dcs->Reset(); | |
efe6147d | 876 | SetBinSizeR(origBinSizeR); |
877 | SetBinSizeZ(origBinSizeZ); | |
6b4d9537 | 878 | return retcode; |
879 | } | |
70c95f95 | 880 | |
f5f6da22 | 881 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
882 | if(validate[i]<1)fLines.RemoveAt(i); | |
70c95f95 | 883 | } |
f5f6da22 | 884 | fLines.Compress(); |
885 | AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast())); | |
70c95f95 | 886 | delete [] validate; |
887 | ||
87104b06 | 888 | // Exit here if Pileup Algorithm 2 has been chosen during second loop |
889 | if(fPileupAlgo == 2 && optCuts==1){ | |
1d60ed44 | 890 | fH3d->Reset(); |
891 | fH3dcs->Reset(); | |
efe6147d | 892 | SetBinSizeR(origBinSizeR); |
893 | SetBinSizeZ(origBinSizeZ); | |
87104b06 | 894 | return 0; |
895 | } | |
6b4d9537 | 896 | |
7340864d | 897 | // Find peaks in histos |
05d1294c | 898 | |
70c95f95 | 899 | Double_t peak[3]={0.,0.,0.}; |
7340864d | 900 | Int_t ntrkl,ntimes; |
1d60ed44 | 901 | FindPeaks(fH3d,peak,ntrkl,ntimes); |
902 | fH3d->Reset(); | |
0496a544 | 903 | Double_t binsizer=(rh-rl)/nbr; |
904 | Double_t binsizez=(zh-zl)/nbz; | |
7d766abb | 905 | if(optCuts==0 && (ntrkl<=2 || ntimes>1)){ |
7340864d | 906 | ntrkl=0; |
907 | ntimes=0; | |
1d60ed44 | 908 | FindPeaks(fH3dcs,peak,ntrkl,ntimes); |
7340864d | 909 | binsizer=(rh-rl)/nbrcs; |
910 | binsizez=(zh-zl)/nbzcs; | |
efe6147d | 911 | if(ntrkl==1 || ntimes>1){ |
1d60ed44 | 912 | fH3dcs->Reset(); |
efe6147d | 913 | SetBinSizeR(origBinSizeR); |
914 | SetBinSizeZ(origBinSizeZ); | |
915 | return retcode; | |
916 | } | |
7340864d | 917 | } |
1d60ed44 | 918 | fH3dcs->Reset(); |
7340864d | 919 | |
efe6147d | 920 | Double_t bs=(binsizer+binsizez)/2.; |
921 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ | |
922 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); | |
923 | if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i); | |
924 | } | |
925 | fLines.Compress(); | |
926 | AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast())); | |
927 | ||
507265fd | 928 | // Finer Histo in limited range in case of high mult. |
efe6147d | 929 | if(rebinned){ |
507265fd | 930 | SetBinSizeR(0.01); |
931 | SetBinSizeZ(0.01); | |
932 | Double_t xl=peak[0]-0.3; | |
933 | Double_t xh=peak[0]+0.3; | |
934 | Double_t yl=peak[1]-0.3; | |
935 | Double_t yh=peak[1]+0.3; | |
936 | zl=peak[2]-0.5; | |
937 | zh=peak[2]+0.5; | |
938 | Int_t nbxfs=(Int_t)((xh-xl)/fBinSizeR+0.0001); | |
939 | Int_t nbyfs=(Int_t)((yh-yl)/fBinSizeR+0.0001); | |
940 | Int_t nbzfs=(Int_t)((zh-zl)/fBinSizeZ+0.0001); | |
941 | ||
1d60ed44 | 942 | if(!fH3dfs){ |
943 | fH3dfs = new TH3F("fH3dfs","xyz distribution",nbxfs,xl,xh,nbyfs,yl,yh,nbzfs,zl,zh); | |
944 | fH3dfs->SetDirectory(0); | |
945 | }else{ | |
946 | fH3dfs->SetBins(nbxfs,xl,xh,nbyfs,yl,yh,nbzfs,zl,zh); | |
947 | } | |
507265fd | 948 | for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){ |
949 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); | |
950 | for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){ | |
951 | AliStrLine *l2 = (AliStrLine*)fLines.At(j); | |
efe6147d | 952 | Double_t dca=l1->GetDCA(l2); |
953 | if(dca > fDCAcut || dca<0.00001) continue; | |
507265fd | 954 | Double_t point[3]; |
955 | Int_t retc = l1->Cross(l2,point); | |
956 | if(retc<0)continue; | |
efe6147d | 957 | Double_t deltaZ=point[2]-zvert; |
958 | if(TMath::Abs(deltaZ)>dZmax)continue; | |
959 | Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]); | |
960 | if(rad>fCoarseMaxRCut)continue; | |
961 | Double_t deltaX=point[0]-xbeam; | |
962 | Double_t deltaY=point[1]-ybeam; | |
963 | Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY); | |
964 | if(raddist>deltaR)continue; | |
1d60ed44 | 965 | fH3dfs->Fill(point[0],point[1],point[2]); |
507265fd | 966 | } |
967 | } | |
968 | ntrkl=0; | |
969 | ntimes=0; | |
970 | ||
971 | Double_t newpeak[3]={0.,0.,0.}; | |
1d60ed44 | 972 | FindPeaks(fH3dfs,newpeak,ntrkl,ntimes); |
507265fd | 973 | if(ntimes==1){ |
974 | for(Int_t iCoo=0; iCoo<3; iCoo++) peak[iCoo]=newpeak[iCoo]; | |
975 | binsizer=fBinSizeR; | |
976 | binsizez=fBinSizeZ; | |
977 | } | |
1d60ed44 | 978 | fH3dfs->Reset(); |
efe6147d | 979 | bs=(binsizer+binsizez)/2.; |
980 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ | |
981 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); | |
982 | if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i); | |
983 | } | |
984 | fLines.Compress(); | |
985 | AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast())); | |
507265fd | 986 | } |
efe6147d | 987 | SetBinSizeR(origBinSizeR); |
988 | SetBinSizeZ(origBinSizeZ); | |
507265fd | 989 | |
990 | ||
70c95f95 | 991 | // Second selection loop |
7340864d | 992 | |
70c95f95 | 993 | |
f5f6da22 | 994 | if(fLines.GetEntriesFast()>1){ |
2ecf7a6f | 995 | retcode=0; |
27167524 | 996 | // find a first candidate for the primary vertex |
f5f6da22 | 997 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); |
70c95f95 | 998 | // make a further selection on tracklets based on this first candidate |
27167524 | 999 | fVert3D.GetXYZ(peak); |
1000 | AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2])); | |
410cfb90 | 1001 | Int_t *validate2 = new Int_t [fLines.GetEntriesFast()]; |
1002 | for(Int_t i=0; i<fLines.GetEntriesFast();i++) validate2[i]=1; | |
f5f6da22 | 1003 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
410cfb90 | 1004 | if(validate2[i]==0) continue; |
f5f6da22 | 1005 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); |
1006 | if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i); | |
410cfb90 | 1007 | if(optCuts==2){ // temporarily only for pileup |
1008 | for(Int_t j=i+1; j<fLines.GetEntriesFast();j++){ | |
1009 | AliStrLine *l2 = (AliStrLine*)fLines.At(j); | |
1010 | if(l1->GetDCA(l2)<0.00001){ | |
1011 | Int_t delta1=(Int_t)l1->GetIdPoint(0)-(Int_t)l2->GetIdPoint(0); | |
1012 | Int_t delta2=(Int_t)l1->GetIdPoint(1)-(Int_t)l2->GetIdPoint(1); | |
1013 | Int_t deltamod1=(Int_t)l1->GetIdPoint(0)/kMaxCluPerMod | |
1014 | -(Int_t)l2->GetIdPoint(0)/kMaxCluPerMod; | |
1015 | Int_t deltamod2=(Int_t)l1->GetIdPoint(1)/kMaxCluPerMod | |
1016 | -(Int_t)l2->GetIdPoint(1)/kMaxCluPerMod; | |
1017 | // remove tracklets sharing a point | |
1018 | if( (delta1==0 && deltamod2==0) || | |
1019 | (delta2==0 && deltamod1==0) ) validate2[j]=0; | |
1020 | } | |
1021 | } | |
1022 | } | |
70c95f95 | 1023 | } |
410cfb90 | 1024 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
1025 | if(validate2[i]==0) fLines.RemoveAt(i); | |
1026 | } | |
1027 | delete [] validate2; | |
f5f6da22 | 1028 | fLines.Compress(); |
1029 | AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast())); | |
2ecf7a6f | 1030 | if(fLines.GetEntriesFast()>1){// this new tracklet selection is used |
1031 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); | |
1032 | } | |
70c95f95 | 1033 | } |
1034 | return retcode; | |
1035 | } | |
1036 | ||
efe6147d | 1037 | //________________________________________________________ |
1038 | Int_t AliITSVertexer3D::Prepare3DVertexPbPb(){ | |
1039 | // Finds the 3D vertex information in Pb-Pb events using tracklets | |
1040 | AliDebug(1,"High multiplicity event.\n"); | |
1041 | ||
1042 | Int_t nxy=(Int_t)(2.*fCoarseMaxRCut/f3DBinSize); | |
1043 | Double_t xymi= -nxy*f3DBinSize/2.; | |
1044 | Double_t xyma= nxy*f3DBinSize/2.; | |
1045 | Int_t nz=(Int_t)(2.*fZCutDiamond/f3DBinSize); | |
1046 | Double_t zmi=-nz*f3DBinSize/2.; | |
1047 | Double_t zma=nz*f3DBinSize/2.; | |
1048 | Int_t nolines=fLines.GetEntriesFast(); | |
1d60ed44 | 1049 | if(!fH3dv)fH3dv = new TH3F("fH3dv","3d tracklets",nxy,xymi,xyma,nxy,xymi,xyma,nz,zmi,zma); |
efe6147d | 1050 | |
1051 | for(Int_t itra=0; itra<nolines; itra++){ | |
1052 | Double_t wei = GetFraction(itra); | |
1053 | //printf("tracklet %d ) - weight %f \n",itra,wei); | |
1054 | if(wei>1.e-6){ | |
1055 | AliStrLine *str=(AliStrLine*)fLines.At(itra); | |
1056 | Double_t t1,t2; | |
1057 | if(str->GetParamAtRadius(fCoarseMaxRCut,t1,t2)){ | |
1058 | do{ | |
1059 | Double_t punt[3]; | |
1060 | str->ComputePointAtT(t1,punt); | |
1d60ed44 | 1061 | fH3dv->Fill(punt[0],punt[1],punt[2],wei); |
efe6147d | 1062 | t1+=f3DBinSize/3.; |
1063 | } while(t1<t2); | |
1064 | } | |
1065 | } | |
1066 | } | |
1067 | Int_t noftrk,noftim; | |
1d60ed44 | 1068 | FindPeaks(fH3dv,f3DPeak,noftrk,noftim); // arg: histo3d, peak, # of contrib., # of other peak with same magnitude |
efe6147d | 1069 | |
1070 | ||
1071 | // Remove all the tracklets which are not passing near peak | |
1072 | ||
1073 | while(nolines--){ | |
1074 | AliStrLine *str=(AliStrLine*)fLines.At(nolines); | |
1075 | Double_t dist = str->GetDistFromPoint(f3DPeak); | |
1076 | if(dist>(2.*f3DBinSize)) fLines.RemoveAt(nolines); | |
1077 | } | |
1078 | fLines.Compress(); | |
1079 | nolines=fLines.GetEntriesFast(); | |
1080 | ||
1d60ed44 | 1081 | fH3dv->Reset(); |
efe6147d | 1082 | |
1083 | Int_t *validate2 = new Int_t [fLines.GetEntriesFast()]; | |
1084 | for(Int_t i=0; i<fLines.GetEntriesFast();i++) validate2[i]=1; | |
1085 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ | |
1086 | if(validate2[i]==0) continue; | |
1087 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); | |
1088 | if(l1->GetDistFromPoint(f3DPeak)> fDCAcut)fLines.RemoveAt(i); | |
1089 | for(Int_t j=i+1; j<fLines.GetEntriesFast();j++){ | |
1090 | AliStrLine *l2 = (AliStrLine*)fLines.At(j); | |
1091 | if(l1->GetDCA(l2)<0.00001){ | |
1092 | Int_t delta1=(Int_t)l1->GetIdPoint(0)-(Int_t)l2->GetIdPoint(0); | |
1093 | Int_t delta2=(Int_t)l1->GetIdPoint(1)-(Int_t)l2->GetIdPoint(1); | |
1094 | Int_t deltamod1=(Int_t)l1->GetIdPoint(0)/kMaxCluPerMod | |
1095 | -(Int_t)l2->GetIdPoint(0)/kMaxCluPerMod; | |
1096 | Int_t deltamod2=(Int_t)l1->GetIdPoint(1)/kMaxCluPerMod | |
1097 | -(Int_t)l2->GetIdPoint(1)/kMaxCluPerMod; | |
1098 | // remove tracklets sharing a point | |
1099 | if( (delta1==0 && deltamod2==0) || | |
1100 | (delta2==0 && deltamod1==0) ) validate2[j]=0; | |
1101 | ||
1102 | } | |
1103 | } | |
1104 | } | |
1105 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ | |
1106 | if(validate2[i]==0) fLines.RemoveAt(i); | |
1107 | } | |
1108 | ||
1109 | delete [] validate2; | |
1110 | fLines.Compress(); | |
1111 | ||
1112 | ||
1113 | AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast())); | |
1114 | ||
1115 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); | |
1116 | fVert3D.GetXYZ(f3DPeak); | |
1117 | ||
1118 | return 0; | |
1119 | } | |
1120 | ||
6a0d56b8 | 1121 | //________________________________________________________ |
0496a544 | 1122 | void AliITSVertexer3D::SetMeanPPtSelTracks(Double_t fieldTesla){ |
5951029f | 1123 | // Sets mean values of Pt based on the field |
1124 | // for P (used in multiple scattering) the most probable value is used | |
6a0d56b8 | 1125 | if(TMath::Abs(fieldTesla-0.5)<0.01){ |
5951029f | 1126 | SetMeanPSelTracks(0.375); |
6a0d56b8 | 1127 | SetMeanPtSelTracks(0.630); |
1128 | }else if(TMath::Abs(fieldTesla-0.4)<0.01){ | |
5951029f | 1129 | SetMeanPSelTracks(0.375); |
6a0d56b8 | 1130 | SetMeanPtSelTracks(0.580); |
1131 | }else if(TMath::Abs(fieldTesla-0.2)<0.01){ | |
5951029f | 1132 | SetMeanPSelTracks(0.375); |
6a0d56b8 | 1133 | SetMeanPtSelTracks(0.530); |
1134 | }else if(fieldTesla<0.00001){ | |
5951029f | 1135 | SetMeanPSelTracks(0.375); |
1136 | SetMeanPtSelTracks(0.230); | |
6a0d56b8 | 1137 | }else{ |
1138 | SetMeanPSelTracks(); | |
1139 | SetMeanPtSelTracks(); | |
70c95f95 | 1140 | } |
70c95f95 | 1141 | } |
1142 | ||
4f40f207 | 1143 | //________________________________________________________ |
1144 | void AliITSVertexer3D::SetZCutDiamond(Double_t zcut){ | |
1145 | // The fiducial region along Z is set. The TH1 object pointed by | |
1146 | // fZHisto is created | |
1147 | if(TMath::AreEqualAbs(zcut,fZCutDiamond,1.e-10))return; | |
1148 | fZCutDiamond=zcut; | |
1149 | if(fZHisto) delete fZHisto; | |
1150 | Double_t binsize=0.02; // default 200 micron | |
1151 | Int_t nbins=static_cast<Int_t>(1+2*fZCutDiamond/binsize); | |
1152 | fZHisto=new TH1F("hz","",nbins,-fZCutDiamond,-fZCutDiamond+binsize*nbins); | |
1153 | fZHisto->SetDirectory(0); | |
1154 | } | |
1155 | ||
7340864d | 1156 | //________________________________________________________ |
1157 | void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){ | |
1158 | // Finds bin with max contents in 3D histo of tracket intersections | |
1159 | TAxis *xax = histo->GetXaxis(); | |
1160 | TAxis *yax = histo->GetYaxis(); | |
1161 | TAxis *zax = histo->GetZaxis(); | |
1162 | peak[0]=0.; | |
1163 | peak[1]=0.; | |
1164 | peak[2]=0.; | |
1165 | nOfTracklets = 0; | |
1166 | nOfTimes=0; | |
7d766abb | 1167 | Int_t peakbin[3]={0,0,0}; |
1168 | Int_t peak2bin[3]={-1,-1,-1}; | |
1169 | Int_t bc2=-1; | |
7340864d | 1170 | for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){ |
242c8025 | 1171 | Double_t xval = xax->GetBinCenter(i); |
7340864d | 1172 | for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){ |
242c8025 | 1173 | Double_t yval = yax->GetBinCenter(j); |
7340864d | 1174 | for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){ |
242c8025 | 1175 | Double_t zval = zax->GetBinCenter(k); |
7340864d | 1176 | Int_t bc =(Int_t)histo->GetBinContent(i,j,k); |
7d766abb | 1177 | if(bc==0) continue; |
7340864d | 1178 | if(bc>nOfTracklets){ |
7d766abb | 1179 | nOfTracklets=bc; |
7340864d | 1180 | peak[2] = zval; |
1181 | peak[1] = yval; | |
1182 | peak[0] = xval; | |
7d766abb | 1183 | peakbin[2] = k; |
1184 | peakbin[1] = j; | |
1185 | peakbin[0] = i; | |
1186 | peak2bin[2] = -1; | |
1187 | peak2bin[1] = -1; | |
1188 | peak2bin[0] = -1; | |
1189 | bc2=-1; | |
7340864d | 1190 | nOfTimes = 1; |
06a7cbee | 1191 | }else if(bc==nOfTracklets){ |
7d766abb | 1192 | if(TMath::Abs(i-peakbin[0])<=1 && TMath::Abs(j-peakbin[1])<=1 && TMath::Abs(k-peakbin[2])<=1){ |
1193 | peak2bin[2] = k; | |
1194 | peak2bin[1] = j; | |
1195 | peak2bin[0] = i; | |
1196 | bc2=bc; | |
1197 | nOfTimes = 1; | |
1198 | }else{ | |
1199 | nOfTimes++; | |
1200 | } | |
7340864d | 1201 | } |
1202 | } | |
1203 | } | |
1204 | } | |
7d766abb | 1205 | if(peak2bin[0]>=-1 && bc2!=-1){ // two contiguous peak-cells with same contents |
1206 | peak[0]=0.5*(xax->GetBinCenter(peakbin[0])+xax->GetBinCenter(peak2bin[0])); | |
1207 | peak[1]=0.5*(yax->GetBinCenter(peakbin[1])+yax->GetBinCenter(peak2bin[1])); | |
1208 | peak[2]=0.5*(zax->GetBinCenter(peakbin[2])+zax->GetBinCenter(peak2bin[2])); | |
1209 | nOfTracklets+=bc2; | |
1210 | nOfTimes=1; | |
1211 | } | |
7340864d | 1212 | } |
cd1c2af1 | 1213 | //________________________________________________________ |
1214 | void AliITSVertexer3D::MarkUsedClusters(){ | |
1215 | // Mark clusters of tracklets used in vertex claulation | |
1216 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ | |
1217 | AliStrLine *lin = (AliStrLine*)fLines.At(i); | |
1218 | Int_t idClu1=lin->GetIdPoint(0); | |
1219 | Int_t idClu2=lin->GetIdPoint(1); | |
1220 | fUsedCluster.SetBitNumber(idClu1); | |
1221 | fUsedCluster.SetBitNumber(idClu2); | |
1222 | } | |
1223 | } | |
1224 | //________________________________________________________ | |
1225 | Int_t AliITSVertexer3D::RemoveTracklets(){ | |
1226 | // Remove trackelts close to first found vertex | |
e690d4d0 | 1227 | Double_t vert[3]={fVert3D.GetX(),fVert3D.GetY(),fVert3D.GetZ()}; |
cd1c2af1 | 1228 | Int_t nRemoved=0; |
1229 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ | |
1230 | AliStrLine *lin = (AliStrLine*)fLines.At(i); | |
1231 | if(lin->GetDistFromPoint(vert)<fDCAforPileup){ | |
1232 | Int_t idClu1=lin->GetIdPoint(0); | |
1233 | Int_t idClu2=lin->GetIdPoint(1); | |
1234 | fUsedCluster.SetBitNumber(idClu1); | |
1235 | fUsedCluster.SetBitNumber(idClu2); | |
1236 | fLines.RemoveAt(i); | |
1237 | ++nRemoved; | |
1238 | } | |
1239 | } | |
1240 | fLines.Compress(); | |
1241 | return nRemoved; | |
1242 | } | |
1243 | //________________________________________________________ | |
1244 | void AliITSVertexer3D::FindOther3DVertices(TTree *itsClusterTree){ | |
1245 | // pileup identification based on 3D vertexing with not used clusters | |
72756d8d | 1246 | |
1247 | fVertArray = new AliESDVertex[kMaxPileupVertices+1]; | |
1248 | fVertArray[0]=(*fCurrentVertex); | |
1249 | Int_t nFoundVert=1; | |
1250 | for(Int_t iPilV=1; iPilV<=kMaxPileupVertices; iPilV++){ | |
1251 | MarkUsedClusters(); | |
1252 | fLines.Clear("C"); | |
1253 | Int_t nolines = FindTracklets(itsClusterTree,2); | |
cd1c2af1 | 1254 | if(nolines>=2){ |
72756d8d | 1255 | Int_t nr=RemoveTracklets(); |
1256 | nolines-=nr; | |
1257 | if(nolines>=2){ | |
1258 | Int_t rc=Prepare3DVertex(2); | |
1259 | if(rc==0){ | |
1260 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); | |
1261 | if(fVert3D.GetNContributors()>=fMinTrackletsForPilup){ | |
1262 | fIsPileup=kTRUE; | |
1263 | fVertArray[nFoundVert]=fVert3D; | |
1264 | nFoundVert++; | |
1265 | if(nFoundVert==2){ | |
e690d4d0 | 1266 | fZpuv=fVert3D.GetZ(); |
72756d8d | 1267 | fNTrpuv=fVert3D.GetNContributors(); |
1268 | } | |
1269 | } | |
cd1c2af1 | 1270 | } |
1271 | } | |
1272 | } | |
1273 | } | |
72756d8d | 1274 | fNoVertices=nFoundVert; |
cd1c2af1 | 1275 | } |
1276 | //______________________________________________________________________ | |
1277 | void AliITSVertexer3D::PileupFromZ(){ | |
1278 | // Calls the pileup algorithm of ALiITSVertexerZ | |
1279 | Int_t binmin, binmax; | |
1280 | Int_t nPeaks=AliITSVertexerZ::GetPeakRegion(fZHisto,binmin,binmax); | |
1281 | if(nPeaks==2)AliWarning("2 peaks found"); | |
1282 | Int_t firstPeakCont=0; | |
242c8025 | 1283 | Double_t firstPeakPos=0.; |
cd1c2af1 | 1284 | for(Int_t i=binmin-1;i<=binmax+1;i++){ |
6b4d9537 | 1285 | firstPeakCont+=static_cast<Int_t>(fZHisto->GetBinContent(i)); |
cd1c2af1 | 1286 | firstPeakPos+=fZHisto->GetBinContent(i)*fZHisto->GetBinCenter(i); |
1287 | } | |
73315c4a | 1288 | if(firstPeakCont>0){ |
1289 | firstPeakPos/=firstPeakCont; | |
1290 | Int_t ncontr2=0; | |
1291 | if(firstPeakCont>fMinTrackletsForPilup){ | |
1292 | Float_t secPeakPos; | |
1293 | ncontr2=AliITSVertexerZ::FindSecondPeak(fZHisto,binmin,binmax,secPeakPos); | |
1294 | if(ncontr2>=fMinTrackletsForPilup){ | |
1295 | fIsPileup=kTRUE; | |
6b4d9537 | 1296 | fNoVertices=2; |
1297 | AliESDVertex secondVert(secPeakPos,0.1,ncontr2); | |
1298 | fVertArray = new AliESDVertex[2]; | |
1299 | fVertArray[0]=(*fCurrentVertex); | |
1300 | fVertArray[1]=secondVert; | |
73315c4a | 1301 | fZpuv=secPeakPos; |
1302 | fNTrpuv=ncontr2; | |
1303 | } | |
cd1c2af1 | 1304 | } |
1305 | } | |
1306 | } | |
efe6147d | 1307 | |
1308 | //________________________________________________________ | |
1309 | Double_t AliITSVertexer3D::GetFraction(Int_t itr) const { | |
1310 | // this method is used to fill a 3D histogram representing | |
1311 | // the trajectories of the candidate tracklets | |
1312 | // The computed fraction is used as a weight at filling time | |
1313 | AliStrLine *str = (AliStrLine*)fLines.At(itr); | |
1314 | Double_t spigolo=10.; | |
1315 | Double_t cd[3]; | |
1316 | str->GetCd(cd); | |
1317 | Double_t par=0.; | |
1318 | Double_t maxl=TMath::Sqrt(3.)*spigolo; | |
1319 | // intersection with a plane normal to the X axis | |
1320 | if(TMath::AreEqualAbs(cd[0],0.,1.e-9)){ | |
1321 | par=1000000.; | |
1322 | } | |
1323 | else { | |
1324 | par=spigolo/cd[0]; | |
1325 | } | |
1326 | Double_t zc=cd[2]*par; | |
1327 | Double_t yc=cd[1]*par; | |
1328 | if((-spigolo<=yc && yc<=spigolo) && (-spigolo<=zc && zc<=spigolo))return TMath::Abs(par/maxl); | |
1329 | // intersection with a plane normal to the Y axis | |
1330 | if(TMath::AreEqualAbs(cd[1],0.,1.e-9)){ | |
1331 | par=1000000.; | |
1332 | } | |
1333 | else { | |
1334 | par=spigolo/cd[1]; | |
1335 | } | |
1336 | zc=cd[2]*par; | |
1337 | Double_t xc=cd[0]*par; | |
1338 | if((-spigolo<=xc && xc<=spigolo) && (-spigolo<=zc && zc<=spigolo))return TMath::Abs(par/maxl); | |
1339 | // intersection with a plane normal to the Z axis | |
1340 | if(TMath::AreEqualAbs(cd[2],0.,1.e-9)){ | |
1341 | par=1000000.; | |
1342 | } | |
1343 | else { | |
1344 | par=spigolo/cd[2]; | |
1345 | } | |
1346 | yc=cd[1]*par; | |
1347 | xc=cd[0]*par; | |
1348 | if((-spigolo<=xc && xc<=spigolo) && (-spigolo<=yc && yc<=spigolo))return TMath::Abs(par/maxl); | |
1349 | // control should never reach the following lines | |
1350 | AliError(Form("anomalous tracklet direction for tracklet %d in fLines\n",itr)); | |
1351 | str->PrintStatus(); | |
1352 | return 0.; | |
1353 | } | |
1354 | ||
70c95f95 | 1355 | //________________________________________________________ |
1356 | void AliITSVertexer3D::PrintStatus() const { | |
1357 | // Print current status | |
87104b06 | 1358 | printf("========= First step selections =====================\n"); |
6a0d56b8 | 1359 | printf("Cut on diamond (Z) %f\n",fZCutDiamond); |
87104b06 | 1360 | printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut); |
1361 | printf("Loose cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut); | |
6a0d56b8 | 1362 | printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut); |
87104b06 | 1363 | printf("========= Second step selections ====================\n"); |
1364 | printf("Cut on tracklet-to-first-vertex Z distance %f\n",fMaxZCut); | |
cd1c2af1 | 1365 | printf("Max Phi difference: %f\n",fDiffPhiMax); |
87104b06 | 1366 | printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut); |
1367 | printf("Cut on tracklet DCA to beam axis (algo2) %f\n",fMaxRCut2); | |
1368 | printf("========= Pileup selections =========================\n"); | |
cd1c2af1 | 1369 | printf("Pileup algo: %d\n",fPileupAlgo); |
87104b06 | 1370 | printf("Min DCA to 1st vertex for pileup (algo 0 and 1): %f\n",fDCAforPileup); |
1371 | printf("Cut on distance between pair-vertices (algo 2): %f\n",fCutOnPairs); | |
efe6147d | 1372 | printf("Maximum number of clusters on L1 or L2 for downscale: %d\n",fMaxNumOfClForDownScale); |
1373 | printf("Maximum number of clusters on L1 or L2 for histo rebin: %d\n",fMaxNumOfClForRebin); | |
6a0d56b8 | 1374 | printf("=======================================================\n"); |
70c95f95 | 1375 | } |
410cfb90 | 1376 |