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