Adding some further mother volumes to speed-up the overlap checking and particle...
[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>
6a0d56b8 16#include "AliESDVertex.h"
cd706e57 17#include "AliLog.h"
6a0d56b8 18#include "AliStrLine.h"
19#include "AliTracker.h"
70c95f95 20#include "AliITSDetTypeRec.h"
21#include "AliITSRecPoint.h"
32e63e47 22#include "AliITSgeomTGeo.h"
6a0d56b8 23#include "AliVertexerTracks.h"
24#include "AliITSVertexer3D.h"
0b7bac55 25#include "AliITSVertexerZ.h"
6b4d9537 26#include "AliITSSortTrkl.h"
70c95f95 27/////////////////////////////////////////////////////////////////
28// this class implements a method to determine
29// the 3 coordinates of the primary vertex
30// for p-p collisions
31// It can be used successfully with Pb-Pb collisions
32////////////////////////////////////////////////////////////////
33
34ClassImp(AliITSVertexer3D)
35
6a0d56b8 36/* $Id$ */
37
70c95f95 38//______________________________________________________________________
39AliITSVertexer3D::AliITSVertexer3D():AliITSVertexer(),
308c2f7c 40fLines("AliStrLine",1000),
70c95f95 41fVert3D(),
42fCoarseDiffPhiCut(0.),
6b4d9537 43fFineDiffPhiCut(0.),
44fCutOnPairs(0.),
05d1294c 45fCoarseMaxRCut(0.),
70c95f95 46fMaxRCut(0.),
47fZCutDiamond(0.),
05d1294c 48fMaxZCut(0.),
70c95f95 49fDCAcut(0.),
6a0d56b8 50fDiffPhiMax(0.),
51fMeanPSelTrk(0.),
cd1c2af1 52fMeanPtSelTrk(0.),
53fUsedCluster(kMaxCluPerMod*kNSPDMod),
54fZHisto(0),
55fDCAforPileup(0.),
56fPileupAlgo(0)
6a0d56b8 57{
70c95f95 58 // Default constructor
59 SetCoarseDiffPhiCut();
6b4d9537 60 SetFineDiffPhiCut();
61 SetCutOnPairs();
05d1294c 62 SetCoarseMaxRCut();
70c95f95 63 SetMaxRCut();
64 SetZCutDiamond();
05d1294c 65 SetMaxZCut();
7203e11a 66 SetDCACut();
70c95f95 67 SetDiffPhiMax();
6a0d56b8 68 SetMeanPSelTracks();
69 SetMeanPtSelTracks();
cd1c2af1 70 SetMinDCAforPileup();
71 SetPileupAlgo();
72 Float_t binsize=0.02; // default 200 micron
73 Int_t nbins=static_cast<Int_t>(1+2*fZCutDiamond/binsize);
74 fZHisto=new TH1F("hz","",nbins,-fZCutDiamond,-fZCutDiamond+binsize*nbins);
70c95f95 75}
76
77//______________________________________________________________________
70c95f95 78AliITSVertexer3D::~AliITSVertexer3D() {
79 // Destructor
f5f6da22 80 fLines.Clear("C");
cd1c2af1 81 if(fZHisto) delete fZHisto;
70c95f95 82}
83
84//______________________________________________________________________
27167524 85void AliITSVertexer3D::ResetVert3D(){
86 //
6b4d9537 87 ResetVertex();
27167524 88 fVert3D.SetXv(0.);
89 fVert3D.SetYv(0.);
90 fVert3D.SetZv(0.);
91 fVert3D.SetDispersion(0.);
92 fVert3D.SetNContributors(0);
cd1c2af1 93 fUsedCluster.ResetAllBits(0);
27167524 94}
95//______________________________________________________________________
308c2f7c 96AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree){
70c95f95 97 // Defines the AliESDVertex for the current event
27167524 98 ResetVert3D();
308c2f7c 99 AliDebug(1,"FindVertexForCurrentEvent - 3D - PROCESSING NEXT EVENT");
6b4d9537 100 fLines.Clear("C");
101 fCurrentVertex = NULL;
70c95f95 102
308c2f7c 103 Int_t nolines = FindTracklets(itsClusterTree,0);
0b7bac55 104 if(nolines>=2){
105 Int_t rc=Prepare3DVertex(0);
6b4d9537 106 if(fPileupAlgo == 2 && rc == 0) FindVertex3DIterative();
107 else if(fPileupAlgo<2 && rc == 0) FindVertex3D(itsClusterTree);
0b7bac55 108 }
6b4d9537 109
0b7bac55 110 if(!fCurrentVertex){
111 AliITSVertexerZ vertz(GetNominalPos()[0],GetNominalPos()[1]);
b96ee725 112 vertz.SetDetTypeRec(GetDetTypeRec());
0b7bac55 113 AliDebug(1,"Call Vertexer Z\n");
cc282e54 114 vertz.SetLowLimit(-fZCutDiamond);
115 vertz.SetHighLimit(fZCutDiamond);
0b7bac55 116 AliESDVertex* vtxz = vertz.FindVertexForCurrentEvent(itsClusterTree);
117 if(vtxz){
118 Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],vtxz->GetZv()};
119 Double_t covmatrix[6];
120 vtxz->GetCovMatrix(covmatrix);
121 Double_t chi2=99999.;
122 Int_t nContr=vtxz->GetNContributors();
123 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
124 fCurrentVertex->SetTitle("vertexer: Z");
125 fCurrentVertex->SetName("SPDVertexZ");
126 delete vtxz;
0b7bac55 127 }
128
70c95f95 129 }
308c2f7c 130 FindMultiplicity(itsClusterTree);
70c95f95 131 return fCurrentVertex;
132}
133
134//______________________________________________________________________
6b4d9537 135void AliITSVertexer3D::FindVertex3D(TTree *itsClusterTree){
136 // 3D algorithm
6b4d9537 137 /* uncomment to debug
138 printf("Vertex found in first iteration:\n");
139 fVert3D.Print();
140 printf("Start second iteration\n");
141 end of debug lines */
142 if(fVert3D.GetNContributors()>0){
143 fLines.Clear("C");
144 Int_t nolines = FindTracklets(itsClusterTree,1);
145 if(nolines>=2){
146 Int_t rc=Prepare3DVertex(1);
2ecf7a6f 147 if(rc!=0) fVert3D.SetNContributors(0); // exclude this vertex
6b4d9537 148 }
149 }
150 /* uncomment to debug
151 printf("Vertex found in second iteration:\n");
152 fVert3D.Print();
153 end of debug lines */
154
155 Float_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv());
156 if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){
157 Double_t position[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
158 Double_t covmatrix[6];
159 fVert3D.GetCovMatrix(covmatrix);
160 Double_t chi2=99999.;
161 Int_t nContr=fVert3D.GetNContributors();
162 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
163 fCurrentVertex->SetTitle("vertexer: 3D");
164 fCurrentVertex->SetName("SPDVertex3D");
165 fCurrentVertex->SetDispersion(fVert3D.GetDispersion());
166 fNoVertices=1;
167
168 switch(fPileupAlgo){
169 case 0: PileupFromZ(); break;
170 case 1: FindOther3DVertices(itsClusterTree); break;
171 default: AliError("Wrong pileup algorithm"); break;
172 }
173 if(fNoVertices==1){
174 fVertArray = new AliESDVertex[1];
175 fVertArray[0]=(*fCurrentVertex);
176 }
177 }
178}
179
180//______________________________________________________________________
181void AliITSVertexer3D::FindVertex3DIterative(){
182 // Defines the AliESDVertex for the current event
183 Int_t numsor=fLines.GetEntriesFast()*(fLines.GetEntriesFast()-1)/2;
184 //cout<<"AliITSVertexer3D::FindVertexForCurentEvent: Number of tracklets selected for vertexing "<<fLines.GetEntriesFast()<<"; Number of pairs: "<<numsor<<endl;
185 AliITSSortTrkl srt(fLines,numsor,fCutOnPairs,fCoarseMaxRCut);
186 srt.FindClusters();
187 AliInfo(Form("Number of vertices: %d",srt.GetNumberOfClusters()));
188
189 fNoVertices = srt.GetNumberOfClusters();
190 //printf("fNoVertices = %d \n",fNoVertices);
191 if(fNoVertices>0){
192 fVertArray = new AliESDVertex[fNoVertices];
193 for(Int_t kk=0; kk<srt.GetNumberOfClusters(); kk++){
194 Int_t size = 0;
195 Int_t *labels = srt.GetTrackletsLab(kk,size);
196 /*
197 Int_t *pairs = srt.GetClusters(kk);
198 Int_t nopai = srt.GetSizeOfCluster(kk);
199 cout<<"***** Vertex number "<<kk<<". Pairs: \n";
200 for(Int_t jj=0;jj<nopai;jj++){
201 cout<<pairs[jj]<<" - ";
202 if(jj>0 & jj%8==0)cout<<endl;
203 }
204 cout<<endl;
205 cout<<"***** Vertex number "<<kk<<". Labels: \n";
206 */
207 AliStrLine **tclo = new AliStrLine* [size];
208 for(Int_t jj=0;jj<size;jj++){
209 // cout<<labels[jj]<<" - ";
210 // if(jj>0 & jj%8==0)cout<<endl;
211 tclo[jj] = dynamic_cast<AliStrLine*>(fLines[labels[jj]]);
212 }
213 // cout<<endl;
214 delete []labels;
215 fVertArray[kk]=AliVertexerTracks::TrackletVertexFinder(tclo,size);
216 delete [] tclo;
217 // fVertArray[kk].PrintStatus();
218 if(kk == 1){
219 // at least one second vertex is present
220 fIsPileup = kTRUE;
221 fNTrpuv = fVertArray[kk].GetNContributors();
222 fZpuv = fVertArray[kk].GetZv();
223 }
224 }
225 Float_t vRadius=TMath::Sqrt(fVertArray[0].GetXv()*fVertArray[0].GetXv()+fVertArray[0].GetYv()*fVertArray[0].GetYv());
226 if(vRadius<GetPipeRadius() && fVertArray[0].GetNContributors()>0){
227 Double_t position[3]={fVertArray[0].GetXv(),fVertArray[0].GetYv(),fVertArray[0].GetZv()};
228 Double_t covmatrix[6];
229 fVertArray[0].GetCovMatrix(covmatrix);
230 Double_t chi2=99999.;
231 Int_t nContr=fVertArray[0].GetNContributors();
232 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr);
233 fCurrentVertex->SetTitle("vertexer: 3D");
234 fCurrentVertex->SetName("SPDVertex3D");
235 fCurrentVertex->SetDispersion(fVertArray[0].GetDispersion());
236 }
237 }
238
239}
240
241//______________________________________________________________________
242Bool_t AliITSVertexer3D::DistBetweenVertices(AliESDVertex &a, AliESDVertex &b, Double_t test, Double_t &dist){
243 // method to compare the distance between vertices a and b with "test"
244 //it returns kTRUE is the distance is less or equal to test
245 dist = (a.GetX()-b.GetX()) * (a.GetX()-b.GetX());
246 dist += (a.GetY()-b.GetY()) * (a.GetY()-b.GetY());
247 dist += (a.GetZ()-b.GetZ()) * (a.GetZ()-b.GetZ());
248 dist = TMath::Sqrt(dist);
249 if(dist <= test)return kTRUE;
250 return kFALSE;
251}
252
253
254//______________________________________________________________________
308c2f7c 255Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){
70c95f95 256 // All the possible combinations between recpoints on layer 1and 2 are
257 // considered. Straight lines (=tracklets)are formed.
258 // The tracklets are processed in Prepare3DVertex
70c95f95 259
6b4d9537 260 if(!GetDetTypeRec())AliFatal("DetTypeRec pointer has not been set");
261
308c2f7c 262 TTree *tR = itsClusterTree;
6b4d9537 263 fDetTypeRec->ResetRecPoints();
06a7cbee 264 fDetTypeRec->SetTreeAddressR(tR);
70c95f95 265 TClonesArray *itsRec = 0;
cd1c2af1 266 if(optCuts==0) fZHisto->Reset();
6b4d9537 267 // gc1 are local and global coordinates for layer 1
6a0d56b8 268 Float_t gc1[3]={0.,0.,0.};
6b4d9537 269 // gc2 are local and global coordinates for layer 2
6a0d56b8 270 Float_t gc2[3]={0.,0.,0.};
70c95f95 271
06a7cbee 272 itsRec = fDetTypeRec->RecPoints();
6b4d9537 273 TBranch *branch = NULL;
70c95f95 274 branch = tR->GetBranch("ITSRecPoints");
6b4d9537 275 if(!branch){
276 AliError("Null pointer for RecPoints branch");
277 return -1;
278 }
70c95f95 279
280 // Set values for cuts
ecce5370 281 Float_t xbeam=GetNominalPos()[0];
282 Float_t ybeam=GetNominalPos()[1];
05d1294c 283 Float_t zvert=0.;
70c95f95 284 Float_t deltaPhi=fCoarseDiffPhiCut;
05d1294c 285 Float_t deltaR=fCoarseMaxRCut;
286 Float_t dZmax=fZCutDiamond;
6b4d9537 287 if(fPileupAlgo == 2){
288 deltaPhi=fFineDiffPhiCut;
289 deltaR=fMaxRCut;
290 if(optCuts != 0)AliWarning(Form("fPileupAlgo=2 AND optCuts=%d has been selected. It should be 0",optCuts));
291 } else if(optCuts==1){
27167524 292 xbeam=fVert3D.GetXv();
293 ybeam=fVert3D.GetYv();
294 zvert=fVert3D.GetZv();
70c95f95 295 deltaPhi = fDiffPhiMax;
05d1294c 296 deltaR=fMaxRCut;
297 dZmax=fMaxZCut;
6b4d9537 298 } else if(optCuts==2){
cd1c2af1 299 xbeam=fVert3D.GetXv();
300 ybeam=fVert3D.GetYv();
301 deltaPhi = fDiffPhiMax;
302 deltaR=fMaxRCut;
70c95f95 303 }
cd1c2af1 304
70c95f95 305 Int_t nrpL1 = 0; // number of rec points on layer 1
306 Int_t nrpL2 = 0; // number of rec points on layer 2
307
308 // By default irstL1=0 and lastL1=79
6a0d56b8 309 Int_t firstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
32e63e47 310 Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;
6a0d56b8 311 for(Int_t module= firstL1; module<=lastL1;module++){ // count number of recopints on layer 1
70c95f95 312 branch->GetEvent(module);
313 nrpL1+= itsRec->GetEntries();
06a7cbee 314 fDetTypeRec->ResetRecPoints();
70c95f95 315 }
6a0d56b8 316 //By default firstL2=80 and lastL2=239
317 Int_t firstL2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
32e63e47 318 Int_t lastL2 = AliITSgeomTGeo::GetModuleIndex(3,1,1)-1;
6a0d56b8 319 for(Int_t module= firstL2; module<=lastL2;module++){ // count number of recopints on layer 2
70c95f95 320 branch->GetEvent(module);
321 nrpL2+= itsRec->GetEntries();
06a7cbee 322 fDetTypeRec->ResetRecPoints();
70c95f95 323 }
324 if(nrpL1 == 0 || nrpL2 == 0){
70c95f95 325 return -1;
326 }
cd706e57 327 AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2));
70c95f95 328
329 Double_t a[3]={xbeam,ybeam,0.};
330 Double_t b[3]={xbeam,ybeam,10.};
331 AliStrLine zeta(a,b,kTRUE);
6a0d56b8 332 Float_t bField=AliTracker::GetBz()/10.; //T
333 SetMeanPPtSelTracks(bField);
70c95f95 334
335 Int_t nolines = 0;
336 // Loop on modules of layer 1
6a0d56b8 337 for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1
8c42830a 338 if(!fUseModule[modul1]) continue;
27167524 339 UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
70c95f95 340 branch->GetEvent(modul1);
341 Int_t nrecp1 = itsRec->GetEntries();
f5f6da22 342 static TClonesArray prpl1("AliITSRecPoint",nrecp1);
343 prpl1.SetOwner();
70c95f95 344 for(Int_t j=0;j<nrecp1;j++){
345 AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
f5f6da22 346 new(prpl1[j])AliITSRecPoint(*recp);
70c95f95 347 }
06a7cbee 348 fDetTypeRec->ResetRecPoints();
70c95f95 349 for(Int_t j=0;j<nrecp1;j++){
cd1c2af1 350 if(j>kMaxCluPerMod) continue;
351 UShort_t idClu1=modul1*kMaxCluPerMod+j;
352 if(fUsedCluster.TestBitNumber(idClu1)) continue;
f5f6da22 353 AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1.At(j);
6a0d56b8 354 recp1->GetGlobalXYZ(gc1);
355 Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam);
70c95f95 356 if(phi1<0)phi1=2*TMath::Pi()+phi1;
27167524 357 for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
358 for(Int_t k=0;k<4;k++){
359 Int_t ladmod=fLadders[ladder-1]+ladl2;
32e63e47 360 if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
361 Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
8c42830a 362 if(!fUseModule[modul2]) continue;
27167524 363 branch->GetEvent(modul2);
364 Int_t nrecp2 = itsRec->GetEntries();
365 for(Int_t j2=0;j2<nrecp2;j2++){
cd1c2af1 366 if(j2>kMaxCluPerMod) continue;
367 UShort_t idClu2=modul2*kMaxCluPerMod+j2;
368 if(fUsedCluster.TestBitNumber(idClu2)) continue;
6a0d56b8 369 AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2);
6a0d56b8 370 recp2->GetGlobalXYZ(gc2);
27167524 371 Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
372 if(phi2<0)phi2=2*TMath::Pi()+phi2;
373 Double_t diff = TMath::Abs(phi2-phi1);
374 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
cd1c2af1 375 if(optCuts==0 && diff<fDiffPhiMax){
376 Float_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
377 Float_t zc1=gc1[2];
378 Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
379 Float_t zc2=gc2[2];
380 Float_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius
381 fZHisto->Fill(zr0);
382 }
27167524 383 if(diff>deltaPhi)continue;
6a0d56b8 384 AliStrLine line(gc1,gc2,kTRUE);
27167524 385 Double_t cp[3];
386 Int_t retcode = line.Cross(&zeta,cp);
387 if(retcode<0)continue;
388 Double_t dca = line.GetDCA(&zeta);
389 if(dca<0.) continue;
390 if(dca>deltaR)continue;
391 Double_t deltaZ=cp[2]-zvert;
392 if(TMath::Abs(deltaZ)>dZmax)continue;
6a0d56b8 393
6a0d56b8 394 if(nolines == 0){
f5f6da22 395 if(fLines.GetEntriesFast()>0)fLines.Clear("C");
6a0d56b8 396 }
397 Float_t cov[6];
398 recp2->GetGlobalCov(cov);
399
cd1c2af1 400
6a0d56b8 401 Float_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]);
402 Float_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
403 Float_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction
404
405 Float_t curvErr=0;
406 if(bField>0.00001){
407 Float_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm
408 Float_t dRad=TMath::Sqrt(TMath::Power((gc1[0]-gc2[0]),2)+TMath::Power((gc1[1]-gc2[1]),2));
409 Float_t aux=dRad/2.+rad1;
410 curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm
411 }
6a0d56b8 412 Float_t sigmasq[3];
413 sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor;
414 sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor;
415 sigmasq[2]=cov[5]*factor*factor;
416
417 // Multiple scattering
418 Float_t beta=1.;
419 Float_t beta2=beta*beta;
420 Float_t p2=fMeanPSelTrk*fMeanPSelTrk;
421 Float_t rBP=GetPipeRadius();
422 Float_t dBP=0.08/35.3; // 800 um of Be
423 Float_t dL1=0.01; //approx. 1% of radiation length
424 Float_t theta2BP=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dBP);
425 Float_t theta2L1=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dL1);
426 Float_t thetaBP=TMath::Sqrt(theta2BP);
427 Float_t thetaL1=TMath::Sqrt(theta2L1);
6a0d56b8 428 for(Int_t ico=0; ico<3;ico++){
429// printf("Error on coord. %d due to cov matrix+curvErr=%f\n",ico,sigmasq[ico]);
430// // sigmasq[ico]+=rad1*rad1*geomfac[ico]*theta2L1/2; // multiple scattering in layer 1
431// // sigmasq[ico]+=rBP*rBP*geomfac[ico]*theta2BP/2; // multiple scattering in beam pipe
432 sigmasq[ico]+=TMath::Power(rad1*TMath::Tan(thetaL1),2)/3.;
433 sigmasq[ico]+=TMath::Power(rBP*TMath::Tan(thetaBP),2)/3.;
434
435// printf("Multipl. scatt. contr %d = %f (LAY1), %f (BP)\n",ico,rad1*rad1*geomfac[ico]*theta2L1/2,rBP*rBP*geomfac[ico]*theta2BP/2);
436// printf("Total error on coord %d = %f\n",ico,sigmasq[ico]);
437 }
438 Float_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.};
439 if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0];
440 if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1];
441 if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2];
cd1c2af1 442 new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE,idClu1,idClu2);
6a0d56b8 443
27167524 444 }
06a7cbee 445 fDetTypeRec->ResetRecPoints();
70c95f95 446 }
70c95f95 447 }
448 }
f5f6da22 449 prpl1.Clear();
70c95f95 450 }
451 if(nolines == 0)return -2;
452 return nolines;
453}
454
455//______________________________________________________________________
05d1294c 456Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
70c95f95 457 // Finds the 3D vertex information using tracklets
458 Int_t retcode = -1;
459
ecce5370 460 Float_t xbeam=GetNominalPos()[0];
461 Float_t ybeam=GetNominalPos()[1];
05d1294c 462 Float_t zvert=0.;
463 Float_t deltaR=fCoarseMaxRCut;
6b4d9537 464 if(fPileupAlgo == 2) {
465 deltaR=fMaxRCut;
466 if(optCuts!=0)AliWarning(Form("fPileupAlgo=2 AND optCuts=%d. It should be 0",optCuts));
467 }
05d1294c 468 Float_t dZmax=fZCutDiamond;
cd1c2af1 469 if(optCuts==1){
27167524 470 xbeam=fVert3D.GetXv();
471 ybeam=fVert3D.GetYv();
472 zvert=fVert3D.GetZv();
05d1294c 473 deltaR=fMaxRCut;
474 dZmax=fMaxZCut;
cd1c2af1 475 }else if(optCuts==2){
476 xbeam=fVert3D.GetXv();
477 ybeam=fVert3D.GetYv();
478 deltaR=fMaxRCut;
05d1294c 479 }
480
70c95f95 481 Int_t nbr=50;
05d1294c 482 Float_t rl=-fCoarseMaxRCut;
483 Float_t rh=fCoarseMaxRCut;
70c95f95 484 Int_t nbz=100;
485 Float_t zl=-fZCutDiamond;
486 Float_t zh=fZCutDiamond;
487 Float_t binsizer=(rh-rl)/nbr;
488 Float_t binsizez=(zh-zl)/nbz;
7340864d 489 Int_t nbrcs=25;
490 Int_t nbzcs=50;
6b4d9537 491 TH3F *h3d = NULL;
492 TH3F *h3dcs = NULL;
493 if(fPileupAlgo !=2){
494 h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
495 h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh);
496 }
70c95f95 497 // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
f5f6da22 498 Int_t *validate = new Int_t [fLines.GetEntriesFast()];
499 for(Int_t i=0; i<fLines.GetEntriesFast();i++)validate[i]=0;
500 for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){
f5f6da22 501 AliStrLine *l1 = (AliStrLine*)fLines.At(i);
502 for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){
503 AliStrLine *l2 = (AliStrLine*)fLines.At(j);
05d1294c 504 Float_t dca=l1->GetDCA(l2);
505 if(dca > fDCAcut || dca<0.00001) continue;
70c95f95 506 Double_t point[3];
507 Int_t retc = l1->Cross(l2,point);
508 if(retc<0)continue;
6a0d56b8 509 Double_t deltaZ=point[2]-zvert;
510 if(TMath::Abs(deltaZ)>dZmax)continue;
70c95f95 511 Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
05d1294c 512 if(rad>fCoarseMaxRCut)continue;
513 Double_t deltaX=point[0]-xbeam;
514 Double_t deltaY=point[1]-ybeam;
05d1294c 515 Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
05d1294c 516 if(raddist>deltaR)continue;
27167524 517 validate[i]=1;
518 validate[j]=1;
6b4d9537 519 if(fPileupAlgo != 2){
520 h3d->Fill(point[0],point[1],point[2]);
521 h3dcs->Fill(point[0],point[1],point[2]);
522 }
70c95f95 523 }
524 }
525
526
527
70c95f95 528 Int_t numbtracklets=0;
f5f6da22 529 for(Int_t i=0; i<fLines.GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++;
6b4d9537 530 if(numbtracklets<2){
531 delete [] validate;
532 if(fPileupAlgo != 2){
533 delete h3d;
534 delete h3dcs;
535 }
536 return retcode;
537 }
70c95f95 538
f5f6da22 539 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
540 if(validate[i]<1)fLines.RemoveAt(i);
70c95f95 541 }
f5f6da22 542 fLines.Compress();
543 AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast()));
70c95f95 544 delete [] validate;
545
6b4d9537 546 // Exit here if Pileup Algorithm 2 has been chosen
547 if(fPileupAlgo == 2)return 0;
548 //
549
550
7340864d 551 // Find peaks in histos
05d1294c 552
70c95f95 553 Double_t peak[3]={0.,0.,0.};
7340864d 554 Int_t ntrkl,ntimes;
555 FindPeaks(h3d,peak,ntrkl,ntimes);
70c95f95 556 delete h3d;
557
7340864d 558 if(optCuts==0 && ntrkl<=2){
559 ntrkl=0;
560 ntimes=0;
561 FindPeaks(h3dcs,peak,ntrkl,ntimes);
562 binsizer=(rh-rl)/nbrcs;
563 binsizez=(zh-zl)/nbzcs;
564 if(ntrkl==1 || ntimes>1){delete h3dcs; return retcode;}
565 }
566 delete h3dcs;
567
70c95f95 568 // Second selection loop
7340864d 569
70c95f95 570 Float_t bs=(binsizer+binsizez)/2.;
f5f6da22 571 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
572 AliStrLine *l1 = (AliStrLine*)fLines.At(i);
573 if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i);
70c95f95 574 }
f5f6da22 575 fLines.Compress();
576 AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast()));
70c95f95 577
f5f6da22 578 if(fLines.GetEntriesFast()>1){
2ecf7a6f 579 retcode=0;
27167524 580 // find a first candidate for the primary vertex
f5f6da22 581 fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
70c95f95 582 // make a further selection on tracklets based on this first candidate
27167524 583 fVert3D.GetXYZ(peak);
584 AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
f5f6da22 585 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
586 AliStrLine *l1 = (AliStrLine*)fLines.At(i);
587 if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i);
70c95f95 588 }
f5f6da22 589 fLines.Compress();
590 AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast()));
2ecf7a6f 591 if(fLines.GetEntriesFast()>1){// this new tracklet selection is used
592 fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
593 }
70c95f95 594 }
595 return retcode;
596}
597
6a0d56b8 598//________________________________________________________
599void AliITSVertexer3D::SetMeanPPtSelTracks(Float_t fieldTesla){
600 // Sets mean values of P and Pt based on the field
601 if(TMath::Abs(fieldTesla-0.5)<0.01){
602 SetMeanPSelTracks(0.885);
603 SetMeanPtSelTracks(0.630);
604 }else if(TMath::Abs(fieldTesla-0.4)<0.01){
605 SetMeanPSelTracks(0.805);
606 SetMeanPtSelTracks(0.580);
607 }else if(TMath::Abs(fieldTesla-0.2)<0.01){
608 SetMeanPSelTracks(0.740);
609 SetMeanPtSelTracks(0.530);
610 }else if(fieldTesla<0.00001){
611 SetMeanPSelTracks(0.730);
612 SetMeanPtSelTracks(0.510);
613 }else{
614 SetMeanPSelTracks();
615 SetMeanPtSelTracks();
70c95f95 616 }
70c95f95 617}
618
7340864d 619//________________________________________________________
620void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){
621 // Finds bin with max contents in 3D histo of tracket intersections
622 TAxis *xax = histo->GetXaxis();
623 TAxis *yax = histo->GetYaxis();
624 TAxis *zax = histo->GetZaxis();
625 peak[0]=0.;
626 peak[1]=0.;
627 peak[2]=0.;
628 nOfTracklets = 0;
629 nOfTimes=0;
630 for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){
631 Float_t xval = xax->GetBinCenter(i);
632 for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){
633 Float_t yval = yax->GetBinCenter(j);
634 for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){
635 Float_t zval = zax->GetBinCenter(k);
636 Int_t bc =(Int_t)histo->GetBinContent(i,j,k);
637 if(bc>nOfTracklets){
638 nOfTracklets = bc;
639 peak[2] = zval;
640 peak[1] = yval;
641 peak[0] = xval;
642 nOfTimes = 1;
06a7cbee 643 }else if(bc==nOfTracklets){
7340864d 644 nOfTimes++;
645 }
646 }
647 }
648 }
7340864d 649}
70c95f95 650
651//________________________________________________________
cd1c2af1 652void AliITSVertexer3D::MarkUsedClusters(){
653 // Mark clusters of tracklets used in vertex claulation
654 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
655 AliStrLine *lin = (AliStrLine*)fLines.At(i);
656 Int_t idClu1=lin->GetIdPoint(0);
657 Int_t idClu2=lin->GetIdPoint(1);
658 fUsedCluster.SetBitNumber(idClu1);
659 fUsedCluster.SetBitNumber(idClu2);
660 }
661}
662//________________________________________________________
663Int_t AliITSVertexer3D::RemoveTracklets(){
664 // Remove trackelts close to first found vertex
665 Double_t vert[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()};
666 Int_t nRemoved=0;
667 for(Int_t i=0; i<fLines.GetEntriesFast();i++){
668 AliStrLine *lin = (AliStrLine*)fLines.At(i);
669 if(lin->GetDistFromPoint(vert)<fDCAforPileup){
670 Int_t idClu1=lin->GetIdPoint(0);
671 Int_t idClu2=lin->GetIdPoint(1);
672 fUsedCluster.SetBitNumber(idClu1);
673 fUsedCluster.SetBitNumber(idClu2);
674 fLines.RemoveAt(i);
675 ++nRemoved;
676 }
677 }
678 fLines.Compress();
679 return nRemoved;
680}
681//________________________________________________________
682void AliITSVertexer3D::FindOther3DVertices(TTree *itsClusterTree){
683 // pileup identification based on 3D vertexing with not used clusters
684 MarkUsedClusters();
685 fLines.Clear("C");
686 Int_t nolines = FindTracklets(itsClusterTree,2);
687 if(nolines>=2){
688 Int_t nr=RemoveTracklets();
689 nolines-=nr;
690 if(nolines>=2){
691 Int_t rc=Prepare3DVertex(2);
692 if(rc==0){
693 fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0);
694 if(fVert3D.GetNContributors()>=fMinTrackletsForPilup){
695 fIsPileup=kTRUE;
6b4d9537 696 fNoVertices=2;
697 fVertArray = new AliESDVertex[2];
698 fVertArray[0]=(*fCurrentVertex);
699 fVertArray[1]=fVert3D;
cd1c2af1 700 fZpuv=fVert3D.GetZv();
701 fNTrpuv=fVert3D.GetNContributors();
702 }
703 }
704 }
705 }
706}
707//______________________________________________________________________
708void AliITSVertexer3D::PileupFromZ(){
709 // Calls the pileup algorithm of ALiITSVertexerZ
710 Int_t binmin, binmax;
711 Int_t nPeaks=AliITSVertexerZ::GetPeakRegion(fZHisto,binmin,binmax);
712 if(nPeaks==2)AliWarning("2 peaks found");
713 Int_t firstPeakCont=0;
714 Float_t firstPeakPos=0.;
715 for(Int_t i=binmin-1;i<=binmax+1;i++){
6b4d9537 716 firstPeakCont+=static_cast<Int_t>(fZHisto->GetBinContent(i));
cd1c2af1 717 firstPeakPos+=fZHisto->GetBinContent(i)*fZHisto->GetBinCenter(i);
718 }
73315c4a 719 if(firstPeakCont>0){
720 firstPeakPos/=firstPeakCont;
721 Int_t ncontr2=0;
722 if(firstPeakCont>fMinTrackletsForPilup){
723 Float_t secPeakPos;
724 ncontr2=AliITSVertexerZ::FindSecondPeak(fZHisto,binmin,binmax,secPeakPos);
725 if(ncontr2>=fMinTrackletsForPilup){
726 fIsPileup=kTRUE;
6b4d9537 727 fNoVertices=2;
728 AliESDVertex secondVert(secPeakPos,0.1,ncontr2);
729 fVertArray = new AliESDVertex[2];
730 fVertArray[0]=(*fCurrentVertex);
731 fVertArray[1]=secondVert;
73315c4a 732 fZpuv=secPeakPos;
733 fNTrpuv=ncontr2;
734 }
cd1c2af1 735 }
736 }
737}
738//________________________________________________________
70c95f95 739void AliITSVertexer3D::PrintStatus() const {
740 // Print current status
6a0d56b8 741 printf("=======================================================\n");
742 printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut);
743 printf("Cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut);
744 printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut);
745 printf("Cut on diamond (Z) %f\n",fZCutDiamond);
746 printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut);
cd1c2af1 747 printf("Max Phi difference: %f\n",fDiffPhiMax);
748 printf("Pileup algo: %d\n",fPileupAlgo);
749 printf("Min DCA to 1st vetrtex for pileup: %f\n",fDCAforPileup);
6a0d56b8 750 printf("=======================================================\n");
70c95f95 751}