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