Adding some further mother volumes to speed-up the overlap checking and particle...
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerCosmics.cxx
CommitLineData
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
52ClassImp(AliITSVertexerCosmics)
53
54//-------------------------------------------------------------------------
55AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
b8ed1a92 56fMaxDistOnOuterLayer(0),
3acc14d5 57fMinDist2Vtxs(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 85AliESDVertex* 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//-------------------------------------------------------------------------
319void 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//-------------------------------------------------------------------------