1 /**************************************************************************
2 * Copyright(c) 2006-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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>
22 #include<TClonesArray.h>
24 #include "AliRunLoader.h"
25 #include "AliITSLoader.h"
26 #include "AliITSDetTypeRec.h"
27 #include "AliITSRecPoint.h"
28 #include "AliITSgeomTGeo.h"
29 /////////////////////////////////////////////////////////////////
30 // this class implements a method to determine
31 // the 3 coordinates of the primary vertex
33 // It can be used successfully with Pb-Pb collisions
34 ////////////////////////////////////////////////////////////////
36 ClassImp(AliITSVertexer3D)
38 //______________________________________________________________________
39 AliITSVertexer3D::AliITSVertexer3D():AliITSVertexer(),
42 fCoarseDiffPhiCut(0.),
50 // Default constructor
51 SetCoarseDiffPhiCut();
60 //______________________________________________________________________
61 AliITSVertexer3D::AliITSVertexer3D(TString fn): AliITSVertexer(fn),
64 fCoarseDiffPhiCut(0.),
72 // Standard constructor
73 fLines = new TClonesArray("AliStrLine",1000);
74 SetCoarseDiffPhiCut();
83 //______________________________________________________________________
84 AliITSVertexer3D::~AliITSVertexer3D() {
92 //______________________________________________________________________
93 void AliITSVertexer3D::ResetVert3D(){
98 fVert3D.SetDispersion(0.);
99 fVert3D.SetNContributors(0);
101 //______________________________________________________________________
102 AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(Int_t evnumber){
103 // Defines the AliESDVertex for the current event
105 AliDebug(1,Form("FindVertexForCurrentEvent - 3D - PROCESSING EVENT %d",evnumber));
106 if(fLines)fLines->Clear();
108 Int_t nolines = FindTracklets(evnumber,0);
110 if(nolines<2)return fCurrentVertex;
111 Int_t rc=Prepare3DVertex(0);
112 if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(fLines,0);
113 /* uncomment to debug
114 printf("Vertex found in first iteration:\n");
116 printf("Start second iteration\n");
117 end of debug lines */
118 if(fVert3D.GetNContributors()>0){
119 if(fLines) fLines->Delete();
120 nolines = FindTracklets(evnumber,1);
122 rc=Prepare3DVertex(1);
123 if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(fLines,0);
126 /* uncomment to debug
127 printf("Vertex found in second iteration:\n");
129 end of debug lines */
131 Float_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv());
132 if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){
133 fCurrentVertex = new AliESDVertex();
134 fCurrentVertex->SetTitle("vertexer: 3D");
135 fCurrentVertex->SetName("Vertex");
136 fCurrentVertex->SetXv(fVert3D.GetXv());
137 fCurrentVertex->SetYv(fVert3D.GetYv());
138 fCurrentVertex->SetZv(fVert3D.GetZv());
139 fCurrentVertex->SetDispersion(fVert3D.GetDispersion());
140 fCurrentVertex->SetNContributors(fVert3D.GetNContributors());
142 FindMultiplicity(evnumber);
143 return fCurrentVertex;
146 //______________________________________________________________________
147 Int_t AliITSVertexer3D::FindTracklets(Int_t evnumber, Int_t optCuts){
148 // All the possible combinations between recpoints on layer 1and 2 are
149 // considered. Straight lines (=tracklets)are formed.
150 // The tracklets are processed in Prepare3DVertex
151 AliRunLoader *rl =AliRunLoader::GetRunLoader();
152 AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
153 // AliITSgeom* geom = itsLoader->GetITSgeom();
154 itsLoader->LoadRecPoints();
155 rl->GetEvent(evnumber);
157 AliITSDetTypeRec detTypeRec;
159 TTree *tR = itsLoader->TreeR();
160 detTypeRec.SetTreeAddressR(tR);
161 TClonesArray *itsRec = 0;
162 // lc and gc are local and global coordinates for layer 1
163 Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
164 Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
165 // lc2 and gc2 are local and global coordinates for layer 2
166 Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
167 Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
169 itsRec = detTypeRec.RecPoints();
171 branch = tR->GetBranch("ITSRecPoints");
173 // Set values for cuts
174 Float_t xbeam=0., ybeam=0.;
176 Float_t deltaPhi=fCoarseDiffPhiCut;
177 Float_t deltaR=fCoarseMaxRCut;
178 Float_t dZmax=fZCutDiamond;
180 xbeam=fVert3D.GetXv();
181 ybeam=fVert3D.GetYv();
182 zvert=fVert3D.GetZv();
183 deltaPhi = fDiffPhiMax;
187 Int_t nrpL1 = 0; // number of rec points on layer 1
188 Int_t nrpL2 = 0; // number of rec points on layer 2
190 // By default irstL1=0 and lastL1=79
191 Int_t irstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
192 Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;
193 for(Int_t module= irstL1; module<=lastL1;module++){ // count number of recopints on layer 1
194 branch->GetEvent(module);
195 nrpL1+= itsRec->GetEntries();
196 detTypeRec.ResetRecPoints();
198 //By default irstL2=80 and lastL2=239
199 Int_t irstL2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
200 Int_t lastL2 = AliITSgeomTGeo::GetModuleIndex(3,1,1)-1;
201 for(Int_t module= irstL2; module<=lastL2;module++){ // count number of recopints on layer 2
202 branch->GetEvent(module);
203 nrpL2+= itsRec->GetEntries();
204 detTypeRec.ResetRecPoints();
206 if(nrpL1 == 0 || nrpL2 == 0){
207 itsLoader->UnloadRecPoints();
210 AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2));
212 Double_t a[3]={xbeam,ybeam,0.};
213 Double_t b[3]={xbeam,ybeam,10.};
214 AliStrLine zeta(a,b,kTRUE);
217 // Loop on modules of layer 1
218 for(Int_t modul1= irstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1
219 UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
220 branch->GetEvent(modul1);
221 Int_t nrecp1 = itsRec->GetEntries();
222 TClonesArray *prpl1 = new TClonesArray("AliITSRecPoint",nrecp1);
224 TClonesArray &rpl1 = *prpl1;
225 for(Int_t j=0;j<nrecp1;j++){
226 AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
227 new(rpl1[j])AliITSRecPoint(*recp);
229 detTypeRec.ResetRecPoints();
230 for(Int_t j=0;j<nrecp1;j++){
231 AliITSRecPoint *recp = (AliITSRecPoint*)prpl1->At(j);
232 // Local coordinates of this recpoint
234 lc[0]=recp->GetDetLocalX();
235 lc[2]=recp->GetDetLocalZ();
236 geom->LtoG(modul1,lc,gc); // global coordinates
238 recp->GetGlobalXYZ(gc);
239 Double_t phi1 = TMath::ATan2(gc[1]-ybeam,gc[0]-xbeam);
240 if(phi1<0)phi1=2*TMath::Pi()+phi1;
241 for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
242 for(Int_t k=0;k<4;k++){
243 Int_t ladmod=fLadders[ladder-1]+ladl2;
244 if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2);
245 Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1);
246 branch->GetEvent(modul2);
247 Int_t nrecp2 = itsRec->GetEntries();
248 for(Int_t j2=0;j2<nrecp2;j2++){
249 recp = (AliITSRecPoint*)itsRec->At(j2);
251 lc2[0]=recp->GetDetLocalX();
252 lc2[2]=recp->GetDetLocalZ();
253 geom->LtoG(modul2,lc2,gc2);
255 recp->GetGlobalXYZ(gc2);
256 Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
257 if(phi2<0)phi2=2*TMath::Pi()+phi2;
258 Double_t diff = TMath::Abs(phi2-phi1);
259 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
260 if(diff>deltaPhi)continue;
261 AliStrLine line(gc,gc2,kTRUE);
263 Int_t retcode = line.Cross(&zeta,cp);
264 if(retcode<0)continue;
265 Double_t dca = line.GetDCA(&zeta);
267 if(dca>deltaR)continue;
268 Double_t deltaZ=cp[2]-zvert;
269 if(TMath::Abs(deltaZ)>dZmax)continue;
270 MakeTracklet(gc,gc2,nolines);
272 detTypeRec.ResetRecPoints();
278 if(nolines == 0)return -2;
282 //______________________________________________________________________
283 void AliITSVertexer3D::FindVertices(){
284 // computes the vertices of the events in the range FirstEvent - LastEvent
285 AliRunLoader *rl = AliRunLoader::GetRunLoader();
286 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
287 itsLoader->ReloadRecPoints();
288 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
290 FindVertexForCurrentEvent(i);
292 WriteCurrentVertex();
298 //______________________________________________________________________
299 Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
300 // Finds the 3D vertex information using tracklets
306 Float_t deltaR=fCoarseMaxRCut;
307 Float_t dZmax=fZCutDiamond;
309 xbeam=fVert3D.GetXv();
310 ybeam=fVert3D.GetYv();
311 zvert=fVert3D.GetZv();
317 Float_t rl=-fCoarseMaxRCut;
318 Float_t rh=fCoarseMaxRCut;
320 Float_t zl=-fZCutDiamond;
321 Float_t zh=fZCutDiamond;
322 Float_t binsizer=(rh-rl)/nbr;
323 Float_t binsizez=(zh-zl)/nbz;
324 TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
326 // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
327 Int_t *validate = new Int_t [fLines->GetEntriesFast()];
328 for(Int_t i=0; i<fLines->GetEntriesFast();i++)validate[i]=0;
329 for(Int_t i=0; i<fLines->GetEntriesFast()-1;i++){
330 if(validate[i]==1)continue;
331 AliStrLine *l1 = (AliStrLine*)fLines->At(i);
332 for(Int_t j=i+1;j<fLines->GetEntriesFast();j++){
333 AliStrLine *l2 = (AliStrLine*)fLines->At(j);
334 Float_t dca=l1->GetDCA(l2);
335 if(dca > fDCAcut || dca<0.00001) continue;
337 Int_t retc = l1->Cross(l2,point);
339 Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
340 if(rad>fCoarseMaxRCut)continue;
341 Double_t deltaX=point[0]-xbeam;
342 Double_t deltaY=point[1]-ybeam;
343 Double_t deltaZ=point[2]-zvert;
344 Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
345 if(TMath::Abs(deltaZ)>dZmax)continue;
346 if(raddist>deltaR)continue;
349 h3d->Fill(point[0],point[1],point[2]);
355 Int_t numbtracklets=0;
356 for(Int_t i=0; i<fLines->GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++;
357 if(numbtracklets<2){delete [] validate; delete h3d; return retcode; }
359 for(Int_t i=0; i<fLines->GetEntriesFast();i++){
360 if(validate[i]<1)fLines->RemoveAt(i);
363 AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines->GetEntriesFast()));
367 // finds peak in histo
368 TAxis *xax = h3d->GetXaxis();
369 TAxis *yax = h3d->GetYaxis();
370 TAxis *zax = h3d->GetZaxis();
371 Double_t peak[3]={0.,0.,0.};
372 Float_t contref = 0.;
373 for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){
374 Float_t xval = xax->GetBinCenter(i);
375 for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){
376 Float_t yval = yax->GetBinCenter(j);
377 for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){
378 Float_t bc = h3d->GetBinContent(i,j,k);
379 Float_t zval = zax->GetBinCenter(k);
391 // Second selection loop
392 Float_t bs=(binsizer+binsizez)/2.;
393 for(Int_t i=0; i<fLines->GetEntriesFast();i++){
394 AliStrLine *l1 = (AliStrLine*)fLines->At(i);
395 if(l1->GetDistFromPoint(peak)>2.5*bs)fLines->RemoveAt(i);
398 AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines->GetEntriesFast()));
400 if(fLines->GetEntriesFast()>1){
401 // find a first candidate for the primary vertex
402 fVert3D=AliVertexerTracks::TrackletVertexFinder(fLines,0);
403 // make a further selection on tracklets based on this first candidate
404 fVert3D.GetXYZ(peak);
405 AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
406 for(Int_t i=0; i<fLines->GetEntriesFast();i++){
407 AliStrLine *l1 = (AliStrLine*)fLines->At(i);
408 if(l1->GetDistFromPoint(peak)> fDCAcut)fLines->RemoveAt(i);
411 AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines->GetEntriesFast()));
412 if(fLines->GetEntriesFast()>1) retcode=0; // this new tracklet selection is used
413 else retcode =1; // the previous tracklet selection will be used
421 //______________________________________________________________________
422 void AliITSVertexer3D::MakeTracklet(Double_t *pA, Double_t *pB, Int_t &nolines) {
424 TClonesArray &lines = *fLines;
426 if(fLines->GetEntriesFast()>0)fLines->Clear();
428 if(fLines->GetEntriesFast()==fLines->GetSize()){
429 Int_t newsize=(Int_t) 1.5*fLines->GetEntriesFast();
430 fLines->Expand(newsize);
433 new(lines[nolines++])AliStrLine(pA,pB,kTRUE);
436 //______________________________________________________________________
437 void AliITSVertexer3D::MakeTracklet(Float_t *pA, Float_t *pB, Int_t &nolines) {// Makes a tracklet
440 for(Int_t i=0;i<3;i++){
444 MakeTracklet(a,b,nolines);
447 //________________________________________________________
448 void AliITSVertexer3D::PrintStatus() const {
449 // Print current status
450 cout <<"=======================================================\n";
451 cout << "Loose cut on Delta Phi "<<fCoarseDiffPhiCut<<endl;
452 cout << "Cut on tracklet DCA to Z axis "<<fCoarseMaxRCut<<endl;
453 cout << "Cut on tracklet DCA to beam axis "<<fMaxRCut<<endl;
454 cout << "Cut on diamond (Z) "<<fZCutDiamond<<endl;
455 cout << "Cut on DCA - tracklet to tracklet and to vertex "<<fDCAcut<<endl;
456 cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
459 cout <<"=======================================================\n";