Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-2007, 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 :
16 : #include <TClonesArray.h>
17 : #include <TTree.h>
18 : #include "AliLog.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"
25 :
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
31 : // in the ITS.
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)
48 : //
49 : // Origin: A.Dainese, andrea.dainese@lnl.infn.it
50 : //-------------------------------------------------------------------------
51 :
52 118 : ClassImp(AliITSVertexerCosmics)
53 :
54 : //-------------------------------------------------------------------------
55 0 : AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
56 0 : fMaxDistOnOuterLayer(0),
57 0 : fMinDist2Vtxs(0)
58 0 : {
59 : // Default constructor
60 0 : SetFirstLastModules(0,0,79);
61 0 : SetFirstLastModules(1,80,239);
62 0 : SetFirstLastModules(2,240,323);
63 0 : SetFirstLastModules(3,324,499);
64 0 : SetFirstLastModules(4,500,1247);
65 0 : SetFirstLastModules(5,1248,2197);
66 : /*
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);
73 : */
74 0 : SetMaxVtxRadius(0,5.5);
75 0 : SetMaxVtxRadius(1,8.5);
76 0 : SetMaxVtxRadius(2,18.5);
77 0 : SetMaxVtxRadius(3,28.5);
78 0 : SetMaxVtxRadius(4,39.5);
79 0 : SetMaxVtxRadius(5,48.5);
80 :
81 0 : SetMaxDistOnOuterLayer();
82 0 : SetMinDist2Vtxs();
83 0 : }
84 : //--------------------------------------------------------------------------
85 : AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(TTree *itsClusterTree)
86 : {
87 : // Defines the AliESDVertex for the current event
88 :
89 0 : fCurrentVertex = 0;
90 :
91 0 : TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
92 0 : itsClusterTree->SetBranchAddress("ITSRecPoints",&recpoints);
93 :
94 0 : Int_t lay,lad,det;
95 :
96 0 : Double_t pos[3]={0.,0.,0.};
97 0 : Double_t err[3]={100.,100.,100.};
98 : Int_t ncontributors = -1;
99 :
100 : // Search for innermost layer with at least two clusters
101 : // on two different modules
102 : Int_t ilayer=0,ilayer2=0;
103 : Int_t nHitModulesSPDinner=0;
104 0 : while(ilayer<AliITSgeomTGeo::GetNLayers()) {
105 0 : if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) {
106 0 : ilayer++;
107 0 : continue;
108 : }
109 : Int_t nHitModules=0;
110 0 : for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer]; imodule++) {
111 0 : itsClusterTree->GetEvent(imodule);
112 0 : AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
113 0 : lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
114 0 : if(lay!=ilayer) AliFatal("Layer mismatch!");
115 0 : if(recpoints->GetEntriesFast()>0) nHitModules++;
116 : }
117 0 : if(ilayer==0) nHitModulesSPDinner=nHitModules;
118 0 : if(nHitModules>=2) break;
119 0 : ilayer++;
120 0 : }
121 :
122 0 : ilayer2=ilayer+1;
123 0 : while(ilayer2<6) {
124 0 : if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer2)) break;
125 0 : ilayer2++;
126 : }
127 :
128 : // try tracklet on SPD2 and point on SPD1
129 0 : if(ilayer==1 && !AliITSReconstructor::GetRecoParam()->GetLayersToSkip(0) &&
130 0 : nHitModulesSPDinner>0) { ilayer=0; ilayer2=1; }
131 :
132 0 : if(ilayer>4 || ilayer2>5) {
133 0 : AliWarning("Not enough clusters");
134 0 : delete recpoints;
135 0 : fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
136 0 : fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
137 0 : fCurrentVertex->SetNContributors(ncontributors);
138 0 : return fCurrentVertex;
139 : }
140 :
141 :
142 : const Int_t arrSize = 1000;
143 0 : Float_t xclInnLay[arrSize],yclInnLay[arrSize],zclInnLay[arrSize],modclInnLay[arrSize];
144 0 : Float_t e2xclInnLay[arrSize],e2yclInnLay[arrSize],e2zclInnLay[arrSize];
145 0 : Float_t e2xclOutLay[arrSize],e2yclOutLay[arrSize],e2zclOutLay[arrSize];
146 : Int_t nclInnLayStored=0;
147 0 : Float_t xclOutLay[arrSize],yclOutLay[arrSize],zclOutLay[arrSize],modclOutLay[arrSize];
148 : Int_t nclOutLayStored=0;
149 : Int_t nRecPoints,nRecPointsInnLay=0;
150 :
151 0 : Float_t gc[3],gcov[6];
152 :
153 : Float_t matchOutLayValue;
154 : Float_t distxyInnLay,distxyInnLayBest=0.;
155 0 : Double_t p1[3],p2[3],p3[3];
156 :
157 : Float_t xvtx,yvtx,zvtx,rvtx;
158 : Int_t i1InnLayBest=-1,i2InnLayBest=-1;
159 :
160 :
161 : // Collect clusters in the selected layer and the outer one
162 0 : for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer2]; imodule++) {
163 0 : itsClusterTree->GetEvent(imodule);
164 0 : AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
165 0 : lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
166 0 : nRecPoints=recpoints->GetEntriesFast();
167 0 : if(imodule<=fLast[ilayer]) nRecPointsInnLay += nRecPoints;
168 : //printf("cosmics: module %d clusters %d\n",imodule,nRecPoints);
169 0 : for(Int_t irp=0; irp<nRecPoints; irp++) {
170 0 : AliITSRecPoint *rp=(AliITSRecPoint*)recpoints->UncheckedAt(irp);
171 : // Local coordinates of this recpoint
172 0 : rp->GetGlobalXYZ(gc);
173 0 : if(lay==ilayer) { // store InnLay clusters
174 0 : xclInnLay[nclInnLayStored]=gc[0];
175 0 : yclInnLay[nclInnLayStored]=gc[1];
176 0 : zclInnLay[nclInnLayStored]=gc[2];
177 0 : rp->GetGlobalCov(gcov);
178 0 : e2xclInnLay[nclInnLayStored]=gcov[0];
179 0 : e2yclInnLay[nclInnLayStored]=gcov[3];
180 0 : e2zclInnLay[nclInnLayStored]=gcov[5];
181 0 : modclInnLay[nclInnLayStored]=imodule;
182 0 : nclInnLayStored++;
183 0 : }
184 0 : if(lay==ilayer2) { // store OutLay clusters
185 0 : xclOutLay[nclOutLayStored]=gc[0];
186 0 : yclOutLay[nclOutLayStored]=gc[1];
187 0 : zclOutLay[nclOutLayStored]=gc[2];
188 0 : rp->GetGlobalCov(gcov);
189 0 : e2xclOutLay[nclOutLayStored]=gcov[0];
190 0 : e2yclOutLay[nclOutLayStored]=gcov[3];
191 0 : e2zclOutLay[nclOutLayStored]=gcov[5];
192 0 : modclOutLay[nclOutLayStored]=imodule;
193 0 : nclOutLayStored++;
194 0 : }
195 0 : if(nclInnLayStored>arrSize || nclOutLayStored>arrSize) {
196 : //AliFatal("More than arrSize clusters per layer");
197 0 : AliWarning("Too many clusters per layer");
198 0 : delete recpoints;
199 0 : fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
200 0 : fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
201 0 : fCurrentVertex->SetNContributors(ncontributors);
202 0 : return fCurrentVertex;
203 : }
204 0 : }// end clusters in a module
205 : }// end modules
206 :
207 :
208 : Int_t i1InnLay,i2InnLay,iOutLay;
209 :
210 : // build fake vertices
211 : //printf("Building tracklets on layer %d\n",ilayer);
212 :
213 : // InnLay - first cluster
214 0 : for(i1InnLay=0; i1InnLay<nclInnLayStored; i1InnLay++) {
215 0 : p1[0]=xclInnLay[i1InnLay];
216 0 : p1[1]=yclInnLay[i1InnLay];
217 0 : p1[2]=zclInnLay[i1InnLay];
218 : // InnLay - second cluster
219 0 : for(i2InnLay=i1InnLay+1; i2InnLay<nclInnLayStored; i2InnLay++) {
220 0 : if(modclInnLay[i1InnLay]==modclInnLay[i2InnLay]) continue;
221 0 : p2[0]=xclInnLay[i2InnLay];
222 0 : p2[1]=yclInnLay[i2InnLay];
223 0 : p2[2]=zclInnLay[i2InnLay];
224 : // look for point on OutLay
225 0 : AliStrLine innLayline(p1,p2,kTRUE);
226 0 : for(iOutLay=0; iOutLay<nclOutLayStored; iOutLay++) {
227 0 : p3[0]=xclOutLay[iOutLay];
228 0 : p3[1]=yclOutLay[iOutLay];
229 0 : p3[2]=zclOutLay[iOutLay];
230 : //printf("(%f,%f) (%f,%f) (%f,%f) %f\n",p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],innLayline.GetDistFromPoint(p3));
231 0 : matchOutLayValue=innLayline.GetDistFromPoint(p3);
232 0 : distxyInnLay = (p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1]);
233 0 : if(matchOutLayValue<fMaxDistOnOuterLayer &&
234 0 : distxyInnLay>distxyInnLayBest) {
235 : //printf("found\n");
236 : distxyInnLayBest=distxyInnLay;
237 : i1InnLayBest=i1InnLay;
238 : i2InnLayBest=i2InnLay;
239 0 : }
240 : }
241 0 : } // InnLay - second cluster
242 : } // InnLay - first cluster
243 :
244 0 : if(i1InnLayBest>-1 && i2InnLayBest>-1) {
245 0 : xvtx = 0.5*(xclInnLay[i1InnLayBest]+xclInnLay[i2InnLayBest]);
246 0 : yvtx = 0.5*(yclInnLay[i1InnLayBest]+yclInnLay[i2InnLayBest]);
247 0 : zvtx = 0.5*(zclInnLay[i1InnLayBest]+zclInnLay[i2InnLayBest]);
248 0 : rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
249 0 : if(rvtx<fMaxVtxRadius[ilayer]) {
250 : ncontributors = ilayer;
251 0 : pos[0] = xvtx;
252 0 : pos[1] = yvtx;
253 0 : pos[2] = zvtx;
254 0 : err[0]=TMath::Sqrt(0.25*(e2xclInnLay[i1InnLayBest]+e2xclInnLay[i2InnLayBest]));
255 0 : err[1]=TMath::Sqrt(0.25*(e2yclInnLay[i1InnLayBest]+e2yclInnLay[i2InnLayBest]));
256 0 : err[2]=TMath::Sqrt(0.25*(e2zclInnLay[i1InnLayBest]+e2zclInnLay[i2InnLayBest]));
257 0 : }
258 :
259 : } else { // give it a try exchanging InnLay and OutLay
260 :
261 : // OutLay - first cluster
262 0 : for(i1InnLay=0; i1InnLay<nclOutLayStored; i1InnLay++) {
263 0 : p1[0]=xclOutLay[i1InnLay];
264 0 : p1[1]=yclOutLay[i1InnLay];
265 0 : p1[2]=zclOutLay[i1InnLay];
266 : // OutLay - second cluster
267 0 : for(i2InnLay=i1InnLay+1; i2InnLay<nclOutLayStored; i2InnLay++) {
268 0 : if(modclOutLay[i1InnLay]==modclOutLay[i2InnLay]) continue;
269 0 : p2[0]=xclOutLay[i2InnLay];
270 0 : p2[1]=yclOutLay[i2InnLay];
271 0 : p2[2]=zclOutLay[i2InnLay];
272 : // look for point on InnLay
273 0 : AliStrLine outLayline(p1,p2,kTRUE);
274 0 : for(iOutLay=0; iOutLay<nclInnLayStored; iOutLay++) {
275 0 : p3[0]=xclInnLay[iOutLay];
276 0 : p3[1]=yclInnLay[iOutLay];
277 0 : p3[2]=zclInnLay[iOutLay];
278 : //printf(" %f\n",InnLayline.GetDistFromPoint(p3));
279 0 : matchOutLayValue=outLayline.GetDistFromPoint(p3);
280 0 : distxyInnLay = (p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1]);
281 0 : if(matchOutLayValue<fMaxDistOnOuterLayer &&
282 0 : distxyInnLay>distxyInnLayBest) {
283 : distxyInnLayBest=distxyInnLay;
284 : i1InnLayBest=i1InnLay;
285 : i2InnLayBest=i2InnLay;
286 0 : }
287 : }
288 0 : } // OutLay - second cluster
289 : } // OutLay - first cluster
290 :
291 0 : if(i1InnLayBest>-1 && i2InnLayBest>-1) {
292 0 : xvtx = 0.5*(xclOutLay[i1InnLayBest]+xclOutLay[i2InnLayBest]);
293 0 : yvtx = 0.5*(yclOutLay[i1InnLayBest]+yclOutLay[i2InnLayBest]);
294 0 : zvtx = 0.5*(zclOutLay[i1InnLayBest]+zclOutLay[i2InnLayBest]);
295 0 : rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
296 0 : if(rvtx<fMaxVtxRadius[ilayer]) {
297 : ncontributors = ilayer2;
298 0 : pos[0] = xvtx;
299 0 : pos[1] = yvtx;
300 0 : pos[2] = zvtx;
301 0 : err[0]=TMath::Sqrt(0.25*(e2xclOutLay[i1InnLayBest]+e2xclOutLay[i2InnLayBest]));
302 0 : err[1]=TMath::Sqrt(0.25*(e2yclOutLay[i1InnLayBest]+e2yclOutLay[i2InnLayBest]));
303 0 : err[2]=TMath::Sqrt(0.25*(e2zclOutLay[i1InnLayBest]+e2zclOutLay[i2InnLayBest]));
304 0 : }
305 : }
306 : } // give it a try exchanging InnLay and OutLay
307 :
308 0 : fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
309 0 : fCurrentVertex->SetTitle("cosmics fake vertex");
310 0 : fCurrentVertex->SetNContributors(ncontributors);
311 : //fCurrentVertex->Print();
312 0 : if(fComputeMultiplicity) FindMultiplicity(itsClusterTree);
313 :
314 0 : delete recpoints;
315 :
316 0 : return fCurrentVertex;
317 0 : }
318 :
319 : //-------------------------------------------------------------------------
320 : void AliITSVertexerCosmics::PrintStatus() const
321 : {
322 : // Print current status
323 0 : printf("=======================================================\n");
324 0 : printf(" fMaxDistOnOuterLayer: %f\n",fMaxDistOnOuterLayer);
325 0 : printf(" fMaxVtxRadius[0]: %f\n",fMaxVtxRadius[0]);
326 0 : printf(" fMinDist2Vtxs: %f\n",fMinDist2Vtxs);
327 0 : printf("=======================================================\n");
328 0 : }
329 : //-------------------------------------------------------------------------
|
