1 /**************************************************************************
2 * Copyright(c) 1998-2007, 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 **************************************************************************/
16 #include <TClonesArray.h>
19 #include "AliESDVertex.h"
20 #include "AliITSgeomTGeo.h"
21 #include "AliITSRecPoint.h"
22 #include "AliITSReconstructor.h"
23 #include "AliITSVertexerCosmics.h"
24 #include "AliStrLine.h"
26 //------------------------------------------------------------------------
27 // This class implements a method to construct a "fake" primary
28 // vertex for cosmic events in which the muon crosses one of 5 inner
29 // ITS layers. A fake primary vertex is needed for the reconstruction,
30 // with e.g. AliITStrackerSA, of the two tracks produced by the muon
32 // We build pairs of clusters on a given layer and define the fake vertex as
33 // the mid-point of the straight line joining the two clusters.
34 // We use the innermost layer that has at least two clusters.
35 // We reject the background by requiring at least one cluster on the outer
36 // layer, closer than fMaxDistOnOuterLayer to the tracklet prolongation.
37 // We can reject (potentially pathological) events with the muon track
38 // tangential to the layer by the requiring the radial position of
39 // the vertex to be smaller than fMaxVtxRadius.
40 // Due to background clusters, more than one vertex per event can
41 // be found. We consider the first found.
42 // The errors on x,y,z of the vertex are calculated as errors on the mean
43 // of clusters coordinates. Non-diag elements of vertex cov. mat. are set to 0.
44 // The number of contributors set in the AliESDVertex object is the
45 // number of the layer on which the tracklet was built; if this number is -1,
46 // the procedure could not find a vertex position and by default
47 // the vertex coordinates are set to (0,0,0) with large errors (100,100,100)
49 // Origin: A.Dainese, andrea.dainese@lnl.infn.it
50 //-------------------------------------------------------------------------
52 ClassImp(AliITSVertexerCosmics)
54 //-------------------------------------------------------------------------
55 AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
56 fMaxDistOnOuterLayer(0),
59 // Default constructor
60 SetFirstLastModules(0,0,79);
61 SetFirstLastModules(1,80,239);
62 SetFirstLastModules(2,240,323);
63 SetFirstLastModules(3,324,499);
64 SetFirstLastModules(4,500,1247);
65 SetFirstLastModules(5,1248,2197);
67 SetMaxVtxRadius(0,3.5);
68 SetMaxVtxRadius(1,6.5);
69 SetMaxVtxRadius(2,14.5);
70 SetMaxVtxRadius(3,23.5);
71 SetMaxVtxRadius(4,37.5);
72 SetMaxVtxRadius(5,42.5);
74 SetMaxVtxRadius(0,5.5);
75 SetMaxVtxRadius(1,8.5);
76 SetMaxVtxRadius(2,18.5);
77 SetMaxVtxRadius(3,28.5);
78 SetMaxVtxRadius(4,39.5);
79 SetMaxVtxRadius(5,48.5);
81 SetMaxDistOnOuterLayer();
84 //--------------------------------------------------------------------------
85 AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(TTree *itsClusterTree)
87 // Defines the AliESDVertex for the current event
91 TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
92 itsClusterTree->SetBranchAddress("ITSRecPoints",&recpoints);
96 Double_t pos[3]={0.,0.,0.};
97 Double_t err[3]={100.,100.,100.};
98 Int_t ncontributors = -1;
100 // Search for innermost layer with at least two clusters
101 // on two different modules
102 Int_t ilayer=0,ilayer2=0;
104 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) {
109 for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer]; imodule++) {
110 itsClusterTree->GetEvent(imodule);
111 AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
112 lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
113 if(lay!=ilayer) AliFatal("Layer mismatch!");
114 if(recpoints->GetEntriesFast()>0) nHitModules++;
116 if(nHitModules>=2) break;
122 if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer2)) break;
126 // try tracklet on SPD2 and point on SPD1
127 if(ilayer==1 && !AliITSReconstructor::GetRecoParam()->GetLayersToSkip(0)) {ilayer=0; ilayer2=1;}
129 if(ilayer>4 || ilayer2>5) {
130 AliWarning("Not enough clusters");
132 fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
133 fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
134 fCurrentVertex->SetNContributors(ncontributors);
135 return fCurrentVertex;
139 const Int_t arrSize = 1000;
140 Float_t xclInnLay[arrSize],yclInnLay[arrSize],zclInnLay[arrSize],modclInnLay[arrSize];
141 Float_t e2xclInnLay[arrSize],e2yclInnLay[arrSize],e2zclInnLay[arrSize];
142 Float_t e2xclOutLay[arrSize],e2yclOutLay[arrSize],e2zclOutLay[arrSize];
143 Int_t nclInnLayStored=0;
144 Float_t xclOutLay[arrSize],yclOutLay[arrSize],zclOutLay[arrSize],modclOutLay[arrSize];
145 Int_t nclOutLayStored=0;
146 Int_t nRecPoints,nRecPointsInnLay=0;
148 Float_t gc[3],gcov[5];
150 Float_t matchOutLayValue;
151 Float_t distxyInnLay,distxyInnLayBest=0.;
152 Double_t p1[3],p2[3],p3[3];
154 Float_t xvtx,yvtx,zvtx,rvtx;
155 Int_t i1InnLayBest=-1,i2InnLayBest=-1;
158 // Collect clusters in the selected layer and the outer one
159 for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer2]; imodule++) {
160 itsClusterTree->GetEvent(imodule);
161 AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
162 lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
163 nRecPoints=recpoints->GetEntriesFast();
164 if(imodule<=fLast[ilayer]) nRecPointsInnLay += nRecPoints;
165 //printf("cosmics: module %d clusters %d\n",imodule,nRecPoints);
166 for(Int_t irp=0; irp<nRecPoints; irp++) {
167 AliITSRecPoint *rp=(AliITSRecPoint*)recpoints->UncheckedAt(irp);
168 // Local coordinates of this recpoint
169 rp->GetGlobalXYZ(gc);
170 if(lay==ilayer) { // store InnLay clusters
171 xclInnLay[nclInnLayStored]=gc[0];
172 yclInnLay[nclInnLayStored]=gc[1];
173 zclInnLay[nclInnLayStored]=gc[2];
174 rp->GetGlobalCov(gcov);
175 e2xclInnLay[nclInnLayStored]=gcov[0];
176 e2yclInnLay[nclInnLayStored]=gcov[3];
177 e2zclInnLay[nclInnLayStored]=gcov[5];
178 modclInnLay[nclInnLayStored]=imodule;
181 if(lay==ilayer2) { // store OutLay clusters
182 xclOutLay[nclOutLayStored]=gc[0];
183 yclOutLay[nclOutLayStored]=gc[1];
184 zclOutLay[nclOutLayStored]=gc[2];
185 rp->GetGlobalCov(gcov);
186 e2xclOutLay[nclOutLayStored]=gcov[0];
187 e2yclOutLay[nclOutLayStored]=gcov[3];
188 e2zclOutLay[nclOutLayStored]=gcov[5];
189 modclOutLay[nclOutLayStored]=imodule;
192 if(nclInnLayStored>arrSize || nclOutLayStored>arrSize) {
193 //AliFatal("More than arrSize clusters per layer");
194 AliWarning("Too many clusters per layer");
196 fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
197 fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
198 fCurrentVertex->SetNContributors(ncontributors);
199 return fCurrentVertex;
201 }// end clusters in a module
205 Int_t i1InnLay,i2InnLay,iOutLay;
207 // build fake vertices
208 printf("Building tracklets on layer %d\n",ilayer);
210 // InnLay - first cluster
211 for(i1InnLay=0; i1InnLay<nclInnLayStored; i1InnLay++) {
212 p1[0]=xclInnLay[i1InnLay];
213 p1[1]=yclInnLay[i1InnLay];
214 p1[2]=zclInnLay[i1InnLay];
215 // InnLay - second cluster
216 for(i2InnLay=i1InnLay+1; i2InnLay<nclInnLayStored; i2InnLay++) {
217 if(modclInnLay[i1InnLay]==modclInnLay[i2InnLay]) continue;
218 p2[0]=xclInnLay[i2InnLay];
219 p2[1]=yclInnLay[i2InnLay];
220 p2[2]=zclInnLay[i2InnLay];
221 // look for point on OutLay
222 AliStrLine innLayline(p1,p2,kTRUE);
223 for(iOutLay=0; iOutLay<nclOutLayStored; iOutLay++) {
224 p3[0]=xclOutLay[iOutLay];
225 p3[1]=yclOutLay[iOutLay];
226 p3[2]=zclOutLay[iOutLay];
227 //printf("(%f,%f) (%f,%f) (%f,%f) %f\n",p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],innLayline.GetDistFromPoint(p3));
228 matchOutLayValue=innLayline.GetDistFromPoint(p3);
229 distxyInnLay = (p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1]);
230 if(matchOutLayValue<fMaxDistOnOuterLayer &&
231 distxyInnLay>distxyInnLayBest) {
233 distxyInnLayBest=distxyInnLay;
234 i1InnLayBest=i1InnLay;
235 i2InnLayBest=i2InnLay;
238 } // InnLay - second cluster
239 } // InnLay - first cluster
241 if(i1InnLayBest>-1 && i2InnLayBest>-1) {
242 xvtx = 0.5*(xclInnLay[i1InnLayBest]+xclInnLay[i2InnLayBest]);
243 yvtx = 0.5*(yclInnLay[i1InnLayBest]+yclInnLay[i2InnLayBest]);
244 zvtx = 0.5*(zclInnLay[i1InnLayBest]+zclInnLay[i2InnLayBest]);
245 rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
246 if(rvtx<fMaxVtxRadius[ilayer]) {
247 ncontributors = ilayer;
251 err[0]=TMath::Sqrt(0.25*(e2xclInnLay[i1InnLayBest]+e2xclInnLay[i2InnLayBest]));
252 err[1]=TMath::Sqrt(0.25*(e2yclInnLay[i1InnLayBest]+e2yclInnLay[i2InnLayBest]));
253 err[2]=TMath::Sqrt(0.25*(e2zclInnLay[i1InnLayBest]+e2zclInnLay[i2InnLayBest]));
256 } else { // give it a try exchanging InnLay and OutLay
258 // OutLay - first cluster
259 for(i1InnLay=0; i1InnLay<nclOutLayStored; i1InnLay++) {
260 p1[0]=xclOutLay[i1InnLay];
261 p1[1]=yclOutLay[i1InnLay];
262 p1[2]=zclOutLay[i1InnLay];
263 // OutLay - second cluster
264 for(i2InnLay=i1InnLay+1; i2InnLay<nclOutLayStored; i2InnLay++) {
265 if(modclOutLay[i1InnLay]==modclOutLay[i2InnLay]) continue;
266 p2[0]=xclOutLay[i2InnLay];
267 p2[1]=yclOutLay[i2InnLay];
268 p2[2]=zclOutLay[i2InnLay];
269 // look for point on InnLay
270 AliStrLine outLayline(p1,p2,kTRUE);
271 for(iOutLay=0; iOutLay<nclInnLayStored; iOutLay++) {
272 p3[0]=xclInnLay[iOutLay];
273 p3[1]=yclInnLay[iOutLay];
274 p3[2]=zclInnLay[iOutLay];
275 //printf(" %f\n",InnLayline.GetDistFromPoint(p3));
276 matchOutLayValue=outLayline.GetDistFromPoint(p3);
277 distxyInnLay = (p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1]);
278 if(matchOutLayValue<fMaxDistOnOuterLayer &&
279 distxyInnLay>distxyInnLayBest) {
280 distxyInnLayBest=distxyInnLay;
281 i1InnLayBest=i1InnLay;
282 i2InnLayBest=i2InnLay;
285 } // OutLay - second cluster
286 } // OutLay - first cluster
288 if(i1InnLayBest>-1 && i2InnLayBest>-1) {
289 xvtx = 0.5*(xclOutLay[i1InnLayBest]+xclOutLay[i2InnLayBest]);
290 yvtx = 0.5*(yclOutLay[i1InnLayBest]+yclOutLay[i2InnLayBest]);
291 zvtx = 0.5*(zclOutLay[i1InnLayBest]+zclOutLay[i2InnLayBest]);
292 rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
293 if(rvtx<fMaxVtxRadius[ilayer]) {
294 ncontributors = ilayer2;
298 err[0]=TMath::Sqrt(0.25*(e2xclOutLay[i1InnLayBest]+e2xclOutLay[i2InnLayBest]));
299 err[1]=TMath::Sqrt(0.25*(e2yclOutLay[i1InnLayBest]+e2yclOutLay[i2InnLayBest]));
300 err[2]=TMath::Sqrt(0.25*(e2zclOutLay[i1InnLayBest]+e2zclOutLay[i2InnLayBest]));
303 } // give it a try exchanging InnLay and OutLay
305 fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
306 fCurrentVertex->SetTitle("cosmics fake vertex");
307 fCurrentVertex->SetNContributors(ncontributors);
308 //fCurrentVertex->Print();
312 return fCurrentVertex;
315 //-------------------------------------------------------------------------
316 void AliITSVertexerCosmics::PrintStatus() const
318 // Print current status
319 printf("=======================================================\n");
320 printf(" fMaxDistOnOuterLayer: %f\n",fMaxDistOnOuterLayer);
321 printf(" fMaxVtxRadius[0]: %f\n",fMaxVtxRadius[0]);
322 printf(" fMinDist2Vtxs: %f\n",fMinDist2Vtxs);
323 printf("=======================================================\n");
325 //-------------------------------------------------------------------------