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 /////////////////////////////////////////////////////////////////
29 // this class implements a method to determine
30 // the 3 coordinates of the primary vertex
32 // It can be used successfully with Pb-Pb collisions
33 ////////////////////////////////////////////////////////////////
35 ClassImp(AliITSVertexer3D)
37 //______________________________________________________________________
38 AliITSVertexer3D::AliITSVertexer3D():AliITSVertexer(),
41 fCoarseDiffPhiCut(0.),
49 // Default constructor
50 SetCoarseDiffPhiCut();
59 //______________________________________________________________________
60 AliITSVertexer3D::AliITSVertexer3D(TString fn): AliITSVertexer(fn),
63 fCoarseDiffPhiCut(0.),
71 // Standard constructor
72 fLines = new TClonesArray("AliStrLine",1000);
73 SetCoarseDiffPhiCut();
82 //______________________________________________________________________
83 AliITSVertexer3D::~AliITSVertexer3D() {
91 //______________________________________________________________________
92 void AliITSVertexer3D::ResetVert3D(){
97 fVert3D.SetDispersion(0.);
98 fVert3D.SetNContributors(0);
100 //______________________________________________________________________
101 AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(Int_t evnumber){
102 // Defines the AliESDVertex for the current event
104 AliDebug(1,Form("FindVertexForCurrentEvent - 3D - PROCESSING EVENT %d",evnumber));
105 if(fLines)fLines->Clear();
107 Int_t nolines = FindTracklets(evnumber,0);
109 if(nolines<2)return fCurrentVertex;
110 Int_t rc=Prepare3DVertex(0);
111 if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(fLines,0);
112 /* uncomment to debug
113 printf("Vertex found in first iteration:\n");
115 printf("Start second iteration\n");
116 end of debug lines */
117 if(fVert3D.GetNContributors()>0){
118 if(fLines) fLines->Delete();
119 nolines = FindTracklets(evnumber,1);
121 rc=Prepare3DVertex(1);
122 if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(fLines,0);
125 /* uncomment to debug
126 printf("Vertex found in second iteration:\n");
128 end of debug lines */
130 Float_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv());
131 if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){
132 fCurrentVertex = new AliESDVertex();
133 fCurrentVertex->SetTitle("vertexer: 3D");
134 fCurrentVertex->SetName("Vertex");
135 fCurrentVertex->SetXv(fVert3D.GetXv());
136 fCurrentVertex->SetYv(fVert3D.GetYv());
137 fCurrentVertex->SetZv(fVert3D.GetZv());
138 fCurrentVertex->SetDispersion(fVert3D.GetDispersion());
139 fCurrentVertex->SetNContributors(fVert3D.GetNContributors());
141 FindMultiplicity(evnumber);
142 return fCurrentVertex;
145 //______________________________________________________________________
146 Int_t AliITSVertexer3D::FindTracklets(Int_t evnumber, Int_t optCuts){
147 // All the possible combinations between recpoints on layer 1and 2 are
148 // considered. Straight lines (=tracklets)are formed.
149 // The tracklets are processed in Prepare3DVertex
150 AliRunLoader *rl =AliRunLoader::GetRunLoader();
151 AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
152 AliITSgeom* geom = itsLoader->GetITSgeom();
153 itsLoader->LoadRecPoints();
154 rl->GetEvent(evnumber);
156 AliITSDetTypeRec detTypeRec;
158 TTree *tR = itsLoader->TreeR();
159 detTypeRec.SetTreeAddressR(tR);
160 TClonesArray *itsRec = 0;
161 // lc and gc are local and global coordinates for layer 1
162 Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
163 Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
164 // lc2 and gc2 are local and global coordinates for layer 2
165 Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
166 Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
168 itsRec = detTypeRec.RecPoints();
170 branch = tR->GetBranch("ITSRecPoints");
172 // Set values for cuts
173 Float_t xbeam=0., ybeam=0.;
175 Float_t deltaPhi=fCoarseDiffPhiCut;
176 Float_t deltaR=fCoarseMaxRCut;
177 Float_t dZmax=fZCutDiamond;
179 xbeam=fVert3D.GetXv();
180 ybeam=fVert3D.GetYv();
181 zvert=fVert3D.GetZv();
182 deltaPhi = fDiffPhiMax;
186 Int_t nrpL1 = 0; // number of rec points on layer 1
187 Int_t nrpL2 = 0; // number of rec points on layer 2
189 // By default irstL1=0 and lastL1=79
190 Int_t irstL1 = geom->GetStartDet(0);
191 Int_t lastL1 = geom->GetModuleIndex(2,1,1)-1;
192 for(Int_t module= irstL1; module<=lastL1;module++){ // count number of recopints on layer 1
193 branch->GetEvent(module);
194 nrpL1+= itsRec->GetEntries();
195 detTypeRec.ResetRecPoints();
197 //By default irstL2=80 and lastL2=239
198 Int_t irstL2 = geom->GetModuleIndex(2,1,1);
199 Int_t lastL2 = geom->GetLastDet(0);
200 for(Int_t module= irstL2; module<=lastL2;module++){ // count number of recopints on layer 2
201 branch->GetEvent(module);
202 nrpL2+= itsRec->GetEntries();
203 detTypeRec.ResetRecPoints();
205 if(nrpL1 == 0 || nrpL2 == 0){
206 itsLoader->UnloadRecPoints();
209 AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2));
211 Double_t a[3]={xbeam,ybeam,0.};
212 Double_t b[3]={xbeam,ybeam,10.};
213 AliStrLine zeta(a,b,kTRUE);
216 // Loop on modules of layer 1
217 for(Int_t modul1= irstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1
218 UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1
219 branch->GetEvent(modul1);
220 Int_t nrecp1 = itsRec->GetEntries();
221 TClonesArray *prpl1 = new TClonesArray("AliITSRecPoint",nrecp1);
223 TClonesArray &rpl1 = *prpl1;
224 for(Int_t j=0;j<nrecp1;j++){
225 AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
226 new(rpl1[j])AliITSRecPoint(*recp);
228 detTypeRec.ResetRecPoints();
229 for(Int_t j=0;j<nrecp1;j++){
230 AliITSRecPoint *recp = (AliITSRecPoint*)prpl1->At(j);
231 // Local coordinates of this recpoint
232 lc[0]=recp->GetDetLocalX();
233 lc[2]=recp->GetDetLocalZ();
234 geom->LtoG(modul1,lc,gc); // global coordinates
235 Double_t phi1 = TMath::ATan2(gc[1]-ybeam,gc[0]-xbeam);
236 if(phi1<0)phi1=2*TMath::Pi()+phi1;
237 for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){
238 for(Int_t k=0;k<4;k++){
239 Int_t ladmod=fLadders[ladder-1]+ladl2;
240 if(ladmod>geom->GetNladders(2)) ladmod=ladmod-geom->GetNladders(2);
241 Int_t modul2=geom->GetModuleIndex(2,ladmod,k+1);
242 branch->GetEvent(modul2);
243 Int_t nrecp2 = itsRec->GetEntries();
244 for(Int_t j2=0;j2<nrecp2;j2++){
245 recp = (AliITSRecPoint*)itsRec->At(j2);
246 lc2[0]=recp->GetDetLocalX();
247 lc2[2]=recp->GetDetLocalZ();
248 geom->LtoG(modul2,lc2,gc2);
249 Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam);
250 if(phi2<0)phi2=2*TMath::Pi()+phi2;
251 Double_t diff = TMath::Abs(phi2-phi1);
252 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
253 if(diff>deltaPhi)continue;
254 AliStrLine line(gc,gc2,kTRUE);
256 Int_t retcode = line.Cross(&zeta,cp);
257 if(retcode<0)continue;
258 Double_t dca = line.GetDCA(&zeta);
260 if(dca>deltaR)continue;
261 Double_t deltaZ=cp[2]-zvert;
262 if(TMath::Abs(deltaZ)>dZmax)continue;
263 MakeTracklet(gc,gc2,nolines);
265 detTypeRec.ResetRecPoints();
271 if(nolines == 0)return -2;
275 //______________________________________________________________________
276 void AliITSVertexer3D::FindVertices(){
277 // computes the vertices of the events in the range FirstEvent - LastEvent
278 AliRunLoader *rl = AliRunLoader::GetRunLoader();
279 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
280 itsLoader->ReloadRecPoints();
281 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
283 FindVertexForCurrentEvent(i);
285 WriteCurrentVertex();
291 //______________________________________________________________________
292 Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){
293 // Finds the 3D vertex information using tracklets
299 Float_t deltaR=fCoarseMaxRCut;
300 Float_t dZmax=fZCutDiamond;
302 xbeam=fVert3D.GetXv();
303 ybeam=fVert3D.GetYv();
304 zvert=fVert3D.GetZv();
310 Float_t rl=-fCoarseMaxRCut;
311 Float_t rh=fCoarseMaxRCut;
313 Float_t zl=-fZCutDiamond;
314 Float_t zh=fZCutDiamond;
315 Float_t binsizer=(rh-rl)/nbr;
316 Float_t binsizez=(zh-zl)/nbz;
317 TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh);
319 // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed)
320 Int_t *validate = new Int_t [fLines->GetEntriesFast()];
321 for(Int_t i=0; i<fLines->GetEntriesFast();i++)validate[i]=0;
322 for(Int_t i=0; i<fLines->GetEntriesFast()-1;i++){
323 if(validate[i]==1)continue;
324 AliStrLine *l1 = (AliStrLine*)fLines->At(i);
325 for(Int_t j=i+1;j<fLines->GetEntriesFast();j++){
326 AliStrLine *l2 = (AliStrLine*)fLines->At(j);
327 Float_t dca=l1->GetDCA(l2);
328 if(dca > fDCAcut || dca<0.00001) continue;
330 Int_t retc = l1->Cross(l2,point);
332 Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
333 if(rad>fCoarseMaxRCut)continue;
334 Double_t deltaX=point[0]-xbeam;
335 Double_t deltaY=point[1]-ybeam;
336 Double_t deltaZ=point[2]-zvert;
337 Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY);
338 if(TMath::Abs(deltaZ)>dZmax)continue;
339 if(raddist>deltaR)continue;
342 h3d->Fill(point[0],point[1],point[2]);
348 Int_t numbtracklets=0;
349 for(Int_t i=0; i<fLines->GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++;
350 if(numbtracklets<2){delete [] validate; delete h3d; return retcode; }
352 for(Int_t i=0; i<fLines->GetEntriesFast();i++){
353 if(validate[i]<1)fLines->RemoveAt(i);
356 AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines->GetEntriesFast()));
360 // finds peak in histo
361 TAxis *xax = h3d->GetXaxis();
362 TAxis *yax = h3d->GetYaxis();
363 TAxis *zax = h3d->GetZaxis();
364 Double_t peak[3]={0.,0.,0.};
365 Float_t contref = 0.;
366 for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){
367 Float_t xval = xax->GetBinCenter(i);
368 for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){
369 Float_t yval = yax->GetBinCenter(j);
370 for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){
371 Float_t bc = h3d->GetBinContent(i,j,k);
372 Float_t zval = zax->GetBinCenter(k);
384 // Second selection loop
385 Float_t bs=(binsizer+binsizez)/2.;
386 for(Int_t i=0; i<fLines->GetEntriesFast();i++){
387 AliStrLine *l1 = (AliStrLine*)fLines->At(i);
388 if(l1->GetDistFromPoint(peak)>2.5*bs)fLines->RemoveAt(i);
391 AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines->GetEntriesFast()));
393 if(fLines->GetEntriesFast()>1){
394 // find a first candidate for the primary vertex
395 fVert3D=AliVertexerTracks::TrackletVertexFinder(fLines,0);
396 // make a further selection on tracklets based on this first candidate
397 fVert3D.GetXYZ(peak);
398 AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2]));
399 for(Int_t i=0; i<fLines->GetEntriesFast();i++){
400 AliStrLine *l1 = (AliStrLine*)fLines->At(i);
401 if(l1->GetDistFromPoint(peak)> fDCAcut)fLines->RemoveAt(i);
404 AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines->GetEntriesFast()));
405 if(fLines->GetEntriesFast()>1) retcode=0; // this new tracklet selection is used
406 else retcode =1; // the previous tracklet selection will be used
414 //______________________________________________________________________
415 void AliITSVertexer3D::MakeTracklet(Double_t *pA, Double_t *pB, Int_t &nolines) {
417 TClonesArray &lines = *fLines;
419 if(fLines->GetEntriesFast()>0)fLines->Clear();
421 if(fLines->GetEntriesFast()==fLines->GetSize()){
422 Int_t newsize=(Int_t) 1.5*fLines->GetEntriesFast();
423 fLines->Expand(newsize);
426 new(lines[nolines++])AliStrLine(pA,pB,kTRUE);
429 //______________________________________________________________________
430 void AliITSVertexer3D::MakeTracklet(Float_t *pA, Float_t *pB, Int_t &nolines) {// Makes a tracklet
433 for(Int_t i=0;i<3;i++){
437 MakeTracklet(a,b,nolines);
440 //________________________________________________________
441 void AliITSVertexer3D::PrintStatus() const {
442 // Print current status
443 cout <<"=======================================================\n";
444 cout << "Loose cut on Delta Phi "<<fCoarseDiffPhiCut<<endl;
445 cout << "Cut on tracklet DCA to Z axis "<<fCoarseMaxRCut<<endl;
446 cout << "Cut on tracklet DCA to beam axis "<<fMaxRCut<<endl;
447 cout << "Cut on diamond (Z) "<<fZCutDiamond<<endl;
448 cout << "Cut on DCA - tracklet to tracklet and to vertex "<<fDCAcut<<endl;
449 cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
452 cout <<"=======================================================\n";