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