Improved definition of cuts for tracklet selection - minor bug fix
[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 **************************************************************************/
15#include <AliESDVertex.h>
16#include <AliITSVertexer3D.h>
17#include <AliStrLine.h>
18#include <AliVertexerTracks.h>
19#include <Riostream.h>
20#include <TH3F.h>
21#include <TTree.h>
22#include<TClonesArray.h>
23#include "AliRunLoader.h"
24#include "AliITSLoader.h"
25#include "AliITSgeom.h"
26#include "AliITSDetTypeRec.h"
27#include "AliITSRecPoint.h"
28/////////////////////////////////////////////////////////////////
29// this class implements a method to determine
30// the 3 coordinates of the primary vertex
31// for p-p collisions
32// It can be used successfully with Pb-Pb collisions
33////////////////////////////////////////////////////////////////
34
35ClassImp(AliITSVertexer3D)
36
37//______________________________________________________________________
38AliITSVertexer3D::AliITSVertexer3D():AliITSVertexer(),
39fLines(),
40fVert3D(),
41fCoarseDiffPhiCut(0.),
05d1294c 42fCoarseMaxRCut(0.),
70c95f95 43fMaxRCut(0.),
44fZCutDiamond(0.),
05d1294c 45fMaxZCut(0.),
70c95f95 46fDCAcut(0.),
47fDiffPhiMax(0.) {
48 // Default constructor
49 SetCoarseDiffPhiCut();
05d1294c 50 SetCoarseMaxRCut();
70c95f95 51 SetMaxRCut();
52 SetZCutDiamond();
05d1294c 53 SetMaxZCut();
70c95f95 54 SetDCAcut();
55 SetDiffPhiMax();
56
57}
58
59//______________________________________________________________________
60AliITSVertexer3D::AliITSVertexer3D(TString fn): AliITSVertexer(fn),
61fLines(),
62fVert3D(),
05d1294c 63fCoarseDiffPhiCut(0.),
64fCoarseMaxRCut(0.),
70c95f95 65fMaxRCut(0.),
66fZCutDiamond(0.),
05d1294c 67fMaxZCut(0.),
70c95f95 68fDCAcut(0.),
69fDiffPhiMax(0.) {
70 // Standard constructor
71 fLines = new TClonesArray("AliStrLine",1000);
72 SetCoarseDiffPhiCut();
05d1294c 73 SetCoarseMaxRCut();
70c95f95 74 SetMaxRCut();
75 SetZCutDiamond();
05d1294c 76 SetMaxZCut();
70c95f95 77 SetDCAcut();
78 SetDiffPhiMax();
79}
80
81//______________________________________________________________________
82AliITSVertexer3D::~AliITSVertexer3D() {
83 // Destructor
84 fVert3D = 0; // this object is not owned by this class
85 if(fLines){
86 fLines->Delete();
87 delete fLines;
88 }
89}
90
91//______________________________________________________________________
92AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(Int_t evnumber){
93 // Defines the AliESDVertex for the current event
94 if (fDebug) cout<<"\n \n FindVertexForCurrentEvent - 3D - PROCESSING EVENT "<<evnumber<<endl;
95 fVert3D = 0;
96 if(fLines)fLines->Clear();
97
98 Int_t nolines = FindTracklets(evnumber,0);
99 fCurrentVertex = 0;
100 if(nolines<2)return fCurrentVertex;
05d1294c 101 Int_t rc=Prepare3DVertex(0);
70c95f95 102 if(rc==0)Find3DVertex();
103 if(fVert3D){
104 if(fLines) fLines->Delete();
105 nolines = FindTracklets(evnumber,1);
d35ad08f 106 if(nolines>=2){
05d1294c 107 rc=Prepare3DVertex(1);
d35ad08f 108 if(rc==0)Find3DVertex();
109 }
70c95f95 110 }
111
112 if(fVert3D){
113 fCurrentVertex = new AliESDVertex();
114 fCurrentVertex->SetTitle("vertexer: 3D");
115 fCurrentVertex->SetName("Vertex");
116 fCurrentVertex->SetXv(fVert3D->GetXv());
117 fCurrentVertex->SetYv(fVert3D->GetYv());
118 fCurrentVertex->SetZv(fVert3D->GetZv());
119 fCurrentVertex->SetDispersion(fVert3D->GetDispersion());
120 fCurrentVertex->SetNContributors(fVert3D->GetNContributors());
121 }
122 FindMultiplicity(evnumber);
123 return fCurrentVertex;
124}
125
126//______________________________________________________________________
127Int_t AliITSVertexer3D::FindTracklets(Int_t evnumber, Int_t optCuts){
128 // All the possible combinations between recpoints on layer 1and 2 are
129 // considered. Straight lines (=tracklets)are formed.
130 // The tracklets are processed in Prepare3DVertex
131 AliRunLoader *rl =AliRunLoader::GetRunLoader();
132 AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
133 AliITSgeom* geom = itsLoader->GetITSgeom();
134 itsLoader->LoadRecPoints();
135 rl->GetEvent(evnumber);
136
137 AliITSDetTypeRec detTypeRec;
138
139 TTree *tR = itsLoader->TreeR();
140 detTypeRec.SetTreeAddressR(tR);
141 TClonesArray *itsRec = 0;
142 // lc and gc are local and global coordinates for layer 1
143 Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
144 Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
145 // lc2 and gc2 are local and global coordinates for layer 2
146 Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
147 Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
148
149 itsRec = detTypeRec.RecPoints();
150 TBranch *branch;
151 branch = tR->GetBranch("ITSRecPoints");
152
153 // Set values for cuts
154 Float_t xbeam=0., ybeam=0.;
05d1294c 155 Float_t zvert=0.;
70c95f95 156 Float_t deltaPhi=fCoarseDiffPhiCut;
05d1294c 157 Float_t deltaR=fCoarseMaxRCut;
158 Float_t dZmax=fZCutDiamond;
70c95f95 159 if(optCuts){
160 xbeam=fVert3D->GetXv();
161 ybeam=fVert3D->GetYv();
05d1294c 162 zvert=fVert3D->GetZv();
70c95f95 163 deltaPhi = fDiffPhiMax;
05d1294c 164 deltaR=fMaxRCut;
165 dZmax=fMaxZCut;
70c95f95 166 }
70c95f95 167 Int_t nrpL1 = 0; // number of rec points on layer 1
168 Int_t nrpL2 = 0; // number of rec points on layer 2
169
170 // By default irstL1=0 and lastL1=79
171 Int_t irstL1 = geom->GetStartDet(0);
172 Int_t lastL1 = geom->GetModuleIndex(2,1,1)-1;
173 for(Int_t module= irstL1; module<=lastL1;module++){ // count number of recopints on layer 1
174 branch->GetEvent(module);
175 nrpL1+= itsRec->GetEntries();
176 detTypeRec.ResetRecPoints();
177 }
178 //By default irstL2=80 and lastL2=239
179 Int_t irstL2 = geom->GetModuleIndex(2,1,1);
180 Int_t lastL2 = geom->GetLastDet(0);
181 for(Int_t module= irstL2; module<=lastL2;module++){ // count number of recopints on layer 2
182 branch->GetEvent(module);
183 nrpL2+= itsRec->GetEntries();
184 detTypeRec.ResetRecPoints();
185 }
186 if(nrpL1 == 0 || nrpL2 == 0){
187 itsLoader->UnloadRecPoints();
188 return -1;
189 }
190 if(fDebug) printf("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2);
191
192 Double_t a[3]={xbeam,ybeam,0.};
193 Double_t b[3]={xbeam,ybeam,10.};
194 AliStrLine zeta(a,b,kTRUE);
195
196 Int_t nolines = 0;
197 // Loop on modules of layer 1
198 for(Int_t modul1= irstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1
199 branch->GetEvent(modul1);
200 Int_t nrecp1 = itsRec->GetEntries();
201 TClonesArray *prpl1 = new TClonesArray("AliITSRecPoint",nrecp1);
202 prpl1->SetOwner();
203 TClonesArray &rpl1 = *prpl1;
204 for(Int_t j=0;j<nrecp1;j++){
205 AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
206 new(rpl1[j])AliITSRecPoint(*recp);
207 }
208 detTypeRec.ResetRecPoints();
209 for(Int_t j=0;j<nrecp1;j++){
210 AliITSRecPoint *recp = (AliITSRecPoint*)prpl1->At(j);
211 // Local coordinates of this recpoint
212 lc[0]=recp->GetDetLocalX();
213 lc[2]=recp->GetDetLocalZ();
214 geom->LtoG(modul1,lc,gc); // global coordinates
215 Double_t phi1 = TMath::ATan2(gc[1]-ybeam,gc[0]-xbeam);
216 if(phi1<0)phi1=2*TMath::Pi()+phi1;
217 for(Int_t modul2= irstL2; modul2<=lastL2;modul2++){
218 branch->GetEvent(modul2);
219 Int_t nrecp2 = itsRec->GetEntries();
220 for(Int_t j2=0;j2<nrecp2;j2++){
221 recp = (AliITSRecPoint*)itsRec->At(j2);
222 lc2[0]=recp->GetDetLocalX();
223 lc2[2]=recp->GetDetLocalZ();
224 geom->LtoG(modul2,lc2,gc2);
225 Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
226 if(phi2<0)phi2=2*TMath::Pi()+phi2;
227 Double_t diff = TMath::Abs(phi2-phi1);
228 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
229 if(diff>deltaPhi)continue;
230 AliStrLine line(gc,gc2,kTRUE);
231 Double_t cp[3];
232 Int_t retcode = line.Cross(&zeta,cp);
233 if(retcode<0)continue;
234 Double_t dca = line.GetDCA(&zeta);
235 if(dca<0.) continue;
05d1294c 236 if(dca>deltaR)continue;
237 Double_t deltaZ=cp[2]-zvert;
238 if(TMath::Abs(deltaZ)>dZmax)continue;
70c95f95 239 MakeTracklet(gc,gc2,nolines);
240 }
241 detTypeRec.ResetRecPoints();
242 }
243 }
244 delete prpl1;
245 }
246 if(nolines == 0)return -2;
247 return nolines;
248}
249
250//______________________________________________________________________
251void AliITSVertexer3D::FindVertices(){
252 // computes the vertices of the events in the range FirstEvent - LastEvent
253 AliRunLoader *rl = AliRunLoader::GetRunLoader();
254 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
255 itsLoader->ReloadRecPoints();
256 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
257 rl->GetEvent(i);
258 FindVertexForCurrentEvent(i);
259 if(fCurrentVertex){
260 WriteCurrentVertex();
261 }
262 }
263}
264
265
266//______________________________________________________________________
267void AliITSVertexer3D::Find3DVertex(){
268 // Determines the vertex position
269 // same algorithm as in AliVertexerTracks::StrLinVertexFinderMinDist(0)
270 // adapted to pure AliStrLine objects
271
272 Int_t knacc = fLines->GetEntriesFast();
273 Double_t (*vectP0)[3]=new Double_t [knacc][3];
274 Double_t (*vectP1)[3]=new Double_t [knacc][3];
275
276 Double_t sum[3][3];
277 Double_t dsum[3]={0,0,0};
278 for(Int_t i=0;i<3;i++)
279 for(Int_t j=0;j<3;j++)sum[i][j]=0;
280 for(Int_t i=0; i<knacc; i++){
281 AliStrLine *line1 = (AliStrLine*)fLines->At(i);
282
283 Double_t p0[3],cd[3];
284 line1->GetP0(p0);
285 line1->GetCd(cd);
286 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
287 vectP0[i][0]=p0[0];
288 vectP0[i][1]=p0[1];
289 vectP0[i][2]=p0[2];
290 vectP1[i][0]=p1[0];
291 vectP1[i][1]=p1[1];
292 vectP1[i][2]=p1[2];
293
294 Double_t matr[3][3];
295 Double_t dknow[3];
296 AliVertexerTracks::GetStrLinDerivMatrix(p0,p1,matr,dknow);
297
298 for(Int_t iii=0;iii<3;iii++){
299 dsum[iii]+=dknow[iii];
300 for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
301 }
302 }
303
304 Double_t vett[3][3];
305 Double_t det=AliVertexerTracks::GetDeterminant3X3(sum);
306 Double_t initPos[3];
307 Double_t sigma = 0.;
308 for(Int_t i=0;i<3;i++)initPos[i]=0.;
309
310 Int_t goodvert=0;
311
312 if(det!=0){
313 for(Int_t zz=0;zz<3;zz++){
314 for(Int_t ww=0;ww<3;ww++){
315 for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
316 }
317 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
318 initPos[zz]=AliVertexerTracks::GetDeterminant3X3(vett)/det;
319 }
320
321 Double_t rvert=TMath::Sqrt(initPos[0]*initPos[0]+initPos[1]*initPos[1]);
322 Double_t zvert=initPos[2];
05d1294c 323 if(rvert<fCoarseMaxRCut && TMath::Abs(zvert)<fZCutDiamond){
70c95f95 324 for(Int_t i=0; i<knacc; i++){
325 Double_t p0[3]={0,0,0},p1[3]={0,0,0};
326 for(Int_t ii=0;ii<3;ii++){
327 p0[ii]=vectP0[i][ii];
328 p1[ii]=vectP1[i][ii];
329 }
330 sigma+=AliVertexerTracks::GetStrLinMinDist(p0,p1,initPos);
331 }
332 sigma=TMath::Sqrt(sigma);
333 goodvert=1;
334 }
335 }
336 delete vectP0;
337 delete vectP1;
338 if(!goodvert){
339 Warning("Find3DVertex","Finder did not succed");
340 for(Int_t i=0;i<3;i++)initPos[i]=0.;
341 sigma=999;
342 knacc=-1;
343 }
344 if(fVert3D){
345 fVert3D-> SetXYZ(initPos);
346 fVert3D->SetDispersion(sigma);
347 fVert3D->SetNContributors(knacc);
348 }else{
349 fVert3D = new AliVertex(initPos,sigma,knacc);
350 }
351 if(fDebug) fVert3D->Print();
352}
353
354
355
356//______________________________________________________________________
05d1294c 357Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
70c95f95 358 // Finds the 3D vertex information using tracklets
359 Int_t retcode = -1;
360
05d1294c 361 Float_t xbeam=0.;
362 Float_t ybeam=0.;
363 Float_t zvert=0.;
364 Float_t deltaR=fCoarseMaxRCut;
365 Float_t dZmax=fZCutDiamond;
366 if(optCuts){
367 xbeam=fVert3D->GetXv();
368 ybeam=fVert3D->GetYv();
369 zvert=fVert3D->GetZv();
370 deltaR=fMaxRCut;
371 dZmax=fMaxZCut;
372 }
373
70c95f95 374 Int_t nbr=50;
05d1294c 375 Float_t rl=-fCoarseMaxRCut;
376 Float_t rh=fCoarseMaxRCut;
70c95f95 377 Int_t nbz=100;
378 Float_t zl=-fZCutDiamond;
379 Float_t zh=fZCutDiamond;
380 Float_t binsizer=(rh-rl)/nbr;
381 Float_t binsizez=(zh-zl)/nbz;
382 TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
383
384
385 // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
386 Int_t *validate = new Int_t [fLines->GetEntriesFast()];
387 for(Int_t i=0; i<fLines->GetEntriesFast();i++)validate[i]=0;
70c95f95 388 for(Int_t i=0; i<fLines->GetEntriesFast()-1;i++){
05d1294c 389 if(validate[i]==1)continue;
70c95f95 390 AliStrLine *l1 = (AliStrLine*)fLines->At(i);
391 for(Int_t j=i+1;j<fLines->GetEntriesFast();j++){
392 AliStrLine *l2 = (AliStrLine*)fLines->At(j);
05d1294c 393 Float_t dca=l1->GetDCA(l2);
394 if(dca > fDCAcut || dca<0.00001) continue;
70c95f95 395 Double_t point[3];
396 Int_t retc = l1->Cross(l2,point);
397 if(retc<0)continue;
398 Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
05d1294c 399 if(rad>fCoarseMaxRCut)continue;
400 Double_t deltaX=point[0]-xbeam;
401 Double_t deltaY=point[1]-ybeam;
402 Double_t deltaZ=point[2]-zvert;
403 Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
404 if(TMath::Abs(deltaZ)>dZmax)continue;
405 if(raddist>deltaR)continue;
70c95f95 406 validate[i]++;
407 validate[j]++;
408 h3d->Fill(point[0],point[1],point[2]);
409 }
410 }
411
412
413
70c95f95 414 Int_t numbtracklets=0;
415 for(Int_t i=0; i<fLines->GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++;
416 if(numbtracklets<2){delete [] validate; delete h3d; return retcode; }
417
418 for(Int_t i=0; i<fLines->GetEntriesFast();i++){
419 if(validate[i]<1)fLines->RemoveAt(i);
420 }
421 fLines->Compress();
422 if (fDebug) cout<<"Number of tracklets (after compress) "<<fLines->GetEntriesFast()<<endl;
423 delete [] validate;
424
05d1294c 425
426 // finds peak in histo
70c95f95 427 TAxis *xax = h3d->GetXaxis();
428 TAxis *yax = h3d->GetYaxis();
429 TAxis *zax = h3d->GetZaxis();
70c95f95 430 Double_t peak[3]={0.,0.,0.};
431 Float_t contref = 0.;
432 for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){
433 Float_t xval = xax->GetBinCenter(i);
434 for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){
435 Float_t yval = yax->GetBinCenter(j);
436 for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){
437 Float_t bc = h3d->GetBinContent(i,j,k);
438 Float_t zval = zax->GetBinCenter(k);
439 if(bc>contref){
440 contref = bc;
441 peak[2] = zval;
442 peak[1] = yval;
443 peak[0] = xval;
444 }
445 }
446 }
447 }
448 delete h3d;
449
450 // Second selection loop
451 Float_t bs=(binsizer+binsizez)/2.;
05d1294c 452 for(Int_t i=0; i<fLines->GetEntriesFast();i++){
70c95f95 453 AliStrLine *l1 = (AliStrLine*)fLines->At(i);
454 if(l1->GetDistFromPoint(peak)>2.5*bs)fLines->RemoveAt(i);
455 }
456 fLines->Compress();
457 if (fDebug) cout<<"Number of tracklets (after 2nd compression) "<<fLines->GetEntriesFast()<<endl;
458
459
460 if(fLines->GetEntriesFast()>1){
461 Find3DVertex(); // find a first candidate for the primary vertex
462 // make a further selection on tracklets based on this first candidate
463 fVert3D->GetXYZ(peak);
05d1294c 464 if (fDebug) cout<<"FIRSTgv Ma V candidate: "<<peak[0]<<"; "<<peak[1]<<"; "<<peak[2]<<endl;
465 for(Int_t i=0; i<fLines->GetEntriesFast();i++){
70c95f95 466 AliStrLine *l1 = (AliStrLine*)fLines->At(i);
467 if(l1->GetDistFromPoint(peak)> fDCAcut)fLines->RemoveAt(i);
468 }
469 fLines->Compress();
470 if (fDebug) cout<<"Number of tracklets (after 3rd compression) "<<fLines->GetEntriesFast()<<endl;
471 if(fLines->GetEntriesFast()>1){
472 retcode = 0; // this last tracklet selection will be used
473 delete fVert3D;
474 fVert3D = 0;
475 }
476 else {
477 retcode =1; // the previous tracklet selection will be used
478 }
479 }
480 else {
481 retcode = 0;
482 }
483 return retcode;
484}
485
486 //______________________________________________________________________
487void AliITSVertexer3D::MakeTracklet(Double_t *pA, Double_t *pB, Int_t &nolines) {
488 // makes a tracklet
489 TClonesArray &lines = *fLines;
490 if(nolines == 0){
491 if(fLines->GetEntriesFast()>0)fLines->Clear();
492 }
493 if(fLines->GetEntriesFast()==fLines->GetSize()){
494 Int_t newsize=(Int_t) 1.5*fLines->GetEntriesFast();
495 fLines->Expand(newsize);
496 }
497
498 new(lines[nolines++])AliStrLine(pA,pB,kTRUE);
499}
500
501 //______________________________________________________________________
502void AliITSVertexer3D::MakeTracklet(Float_t *pA, Float_t *pB, Int_t &nolines) {// Makes a tracklet
503 //
504 Double_t a[3],b[3];
505 for(Int_t i=0;i<3;i++){
506 a[i] = pA[i];
507 b[i] = pB[i];
508 }
509 MakeTracklet(a,b,nolines);
510}
511
512//________________________________________________________
513void AliITSVertexer3D::PrintStatus() const {
514 // Print current status
515 cout <<"=======================================================\n";
516 cout << "Loose cut on Delta Phi "<<fCoarseDiffPhiCut<<endl;
05d1294c 517 cout << "Cut on tracklet DCA to Z axis "<<fCoarseMaxRCut<<endl;
518 cout << "Cut on tracklet DCA to beam axis "<<fMaxRCut<<endl;
70c95f95 519 cout << "Cut on diamond (Z) "<<fZCutDiamond<<endl;
520 cout << "Cut on DCA - tracklet to tracklet and to vertex "<<fDCAcut<<endl;
521 cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
522
523
524 cout <<"=======================================================\n";
525}