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