MCH DCS aliases removed from production configuration
[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"
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);
66 SetMaxVtxRadius(0,3.5);
67 SetMaxVtxRadius(1,6.5);
68 SetMaxVtxRadius(2,14.5);
69 SetMaxVtxRadius(3,23.5);
70 SetMaxVtxRadius(4,37.5);
71 SetMaxVtxRadius(5,42.5);
72 SetMaxDistOnOuterLayer();
3acc14d5 73 SetMinDist2Vtxs();
74}
75//--------------------------------------------------------------------------
76AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(Int_t evnumber)
77{
78 // Defines the AliESDVertex for the current event
79
80 fCurrentVertex = 0;
81 AliRunLoader *rl =AliRunLoader::GetRunLoader();
82 AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
3acc14d5 83 itsLoader->LoadRecPoints();
84 rl->GetEvent(evnumber);
85
86 TTree *rpTree = itsLoader->TreeR();
87
88 TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
89 rpTree->SetBranchAddress("ITSRecPoints",&recpoints);
90
3acc14d5 91 Int_t lay,lad,det;
92
b8ed1a92 93 // Search for innermost layer with at least two clusters
94 // on two different modules
95 Int_t ilayer=0;
96 while(ilayer<6) {
97 Int_t nHitModules=0;
98 for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer]; imodule++) {
99 rpTree->GetEvent(imodule);
100 AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
101 lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
102 if(lay!=ilayer) AliFatal("Layer mismatch!");
103 if(recpoints->GetEntriesFast()>0) nHitModules++;
104 }
105 if(nHitModules>=2) break;
106 ilayer++;
107 }
108 printf("Building tracklets on layer %d\n",ilayer);
109
110
111 Float_t xclInnLay[100],yclInnLay[100],zclInnLay[100],modclInnLay[100];
112 Float_t e2xclInnLay[100],e2yclInnLay[100],e2zclInnLay[100];
113 Int_t nclInnLayStored=0;
114 Float_t xclOutLay[100],yclOutLay[100],zclOutLay[100],modclOutLay[100];
115 Int_t nclOutLayStored=0;
116 Int_t nRecPoints,nRecPointsInnLay=0;
117
118 Float_t gc[3],gcov[5];
119
120 Float_t x[100],y[100],z[100],e2x[100],e2y[100],e2z[100];
121 Double_t p1[3],p2[3],p3[3];
3acc14d5 122 Int_t nvtxs;
b8ed1a92 123 Bool_t good,matchtoOutLay;
124 Float_t xvtx,yvtx,zvtx,rvtx;
3acc14d5 125
b8ed1a92 126 // Collect clusters in the selected layer and the outer one
127 for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer+1]; imodule++) {
3acc14d5 128 rpTree->GetEvent(imodule);
0333e459 129 AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
b8ed1a92 130 lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
131 nRecPoints=recpoints->GetEntriesFast();
132 if(imodule<=fLast[ilayer]) nRecPointsInnLay += nRecPoints;
133 //printf("cosmics: module %d clusters %d\n",imodule,nRecPoints);
134 for(Int_t irp=0; irp<nRecPoints; irp++) {
3acc14d5 135 AliITSRecPoint *rp=(AliITSRecPoint*)recpoints->UncheckedAt(irp);
136 // Local coordinates of this recpoint
b8ed1a92 137 rp->GetGlobalXYZ(gc);
138 if(lay==ilayer) { // store InnLay clusters
139 xclInnLay[nclInnLayStored]=gc[0];
140 yclInnLay[nclInnLayStored]=gc[1];
141 zclInnLay[nclInnLayStored]=gc[2];
142 rp->GetGlobalCov(gcov);
143 e2xclInnLay[nclInnLayStored]=gcov[0];
144 e2yclInnLay[nclInnLayStored]=gcov[3];
145 e2zclInnLay[nclInnLayStored]=gcov[5];
146 modclInnLay[nclInnLayStored]=imodule;
147 nclInnLayStored++;
3acc14d5 148 }
b8ed1a92 149 if(lay==ilayer+1) { // store OutLay clusters
150 xclOutLay[nclOutLayStored]=gc[0];
151 yclOutLay[nclOutLayStored]=gc[1];
152 zclOutLay[nclOutLayStored]=gc[2];
153 modclOutLay[nclOutLayStored]=imodule;
154 nclOutLayStored++;
3acc14d5 155 }
b8ed1a92 156 if(nclInnLayStored>100 || nclOutLayStored>100)
157 AliFatal("More than 100 clusters per layer");
3acc14d5 158 }// end clusters in a module
b8ed1a92 159 }// end modules
3acc14d5 160
161 // build fake vertices
162 nvtxs=0;
b8ed1a92 163 // InnLay - first cluster
164 for(Int_t i1InnLay=0; i1InnLay<nclInnLayStored; i1InnLay++) {
165 p1[0]=xclInnLay[i1InnLay];
166 p1[1]=yclInnLay[i1InnLay];
167 p1[2]=zclInnLay[i1InnLay];
168 // InnLay - second cluster
169 for(Int_t i2InnLay=i1InnLay+1; i2InnLay<nclInnLayStored; i2InnLay++) {
170 if(modclInnLay[i1InnLay]==modclInnLay[i2InnLay]) continue;
171 p2[0]=xclInnLay[i2InnLay];
172 p2[1]=yclInnLay[i2InnLay];
173 p2[2]=zclInnLay[i2InnLay];
174 // look for point on OutLay
175 AliStrLine InnLayline(p1,p2,kTRUE);
176 matchtoOutLay = kFALSE;
177 for(Int_t iOutLay=0; iOutLay<nclOutLayStored; iOutLay++) {
178 p3[0]=xclOutLay[iOutLay];
179 p3[1]=yclOutLay[iOutLay];
180 p3[2]=zclOutLay[iOutLay];
181 //printf(" %f\n",InnLayline.GetDistFromPoint(p3));
182 if(InnLayline.GetDistFromPoint(p3)<fMaxDistOnOuterLayer)
183 { matchtoOutLay = kTRUE; break; }
3acc14d5 184 }
b8ed1a92 185 if(!matchtoOutLay) continue;
186 xvtx = 0.5*(xclInnLay[i1InnLay]+xclInnLay[i2InnLay]);
187 yvtx = 0.5*(yclInnLay[i1InnLay]+yclInnLay[i2InnLay]);
188 zvtx = 0.5*(zclInnLay[i1InnLay]+zclInnLay[i2InnLay]);
3acc14d5 189 rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
b8ed1a92 190 if(rvtx>fMaxVtxRadius[ilayer]) continue;
3acc14d5 191 good = kTRUE;
192 for(Int_t iv=0; iv<nvtxs; iv++) {
193 if(TMath::Sqrt((xvtx- x[iv])*(xvtx- x[iv])+
194 (yvtx- y[iv])*(yvtx- y[iv])+
195 (zvtx- z[iv])*(zvtx- z[iv])) < fMinDist2Vtxs)
196 good = kFALSE;
197 }
198 if(good) {
199 x[nvtxs]=xvtx;
200 y[nvtxs]=yvtx;
201 z[nvtxs]=zvtx;
b8ed1a92 202 e2x[nvtxs]=0.25*(e2xclInnLay[i1InnLay]+e2xclInnLay[i2InnLay]);
203 e2y[nvtxs]=0.25*(e2yclInnLay[i1InnLay]+e2yclInnLay[i2InnLay]);
204 e2z[nvtxs]=0.25*(e2zclInnLay[i1InnLay]+e2zclInnLay[i2InnLay]);
3acc14d5 205 nvtxs++;
206 }
b8ed1a92 207 } // InnLay - second cluster
208 } // InnLay - first cluster
3acc14d5 209
210
211 Double_t pos[3]={0.,0.,0.};
212 Double_t err[3]={100.,100.,100.};
213 if(nvtxs) {
214 pos[0]=x[0];
215 pos[1]=y[0];
216 pos[2]=z[0];
b8ed1a92 217 err[0]=TMath::Sqrt(e2x[0]);
218 err[1]=TMath::Sqrt(e2y[0]);
219 err[2]=TMath::Sqrt(e2z[0]);
3acc14d5 220 }
3acc14d5 221 fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
b8ed1a92 222 if(nvtxs) {
223 fCurrentVertex->SetNContributors(ilayer);
224 } else {
225 fCurrentVertex->SetNContributors(-1);
226 }
3acc14d5 227 fCurrentVertex->SetTitle("cosmics fake vertex");
228
b8ed1a92 229 if(nvtxs>=0) fCurrentVertex->Print();
3acc14d5 230
231 delete recpoints;
232 itsLoader->UnloadRecPoints();
233
234 return fCurrentVertex;
235}
236//-------------------------------------------------------------------------
237void AliITSVertexerCosmics::FindVertices()
238{
239 // computes the vertices of the events in the range FirstEvent - LastEvent
240 AliRunLoader *rl = AliRunLoader::GetRunLoader();
241 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
242 itsLoader->ReloadRecPoints();
243 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
244 // printf("Processing event %d\n",i);
245 rl->GetEvent(i);
246 FindVertexForCurrentEvent(i);
247 if(fCurrentVertex){
248 WriteCurrentVertex();
249 }
250 }
251}
252//-------------------------------------------------------------------------
253void AliITSVertexerCosmics::PrintStatus() const
254{
255 // Print current status
256 printf("=======================================================\n");
b8ed1a92 257 printf(" fMaxDistOnOuterLayer: %f\n",fMaxDistOnOuterLayer);
258 printf(" fMaxVtxRadius[0]: %f\n",fMaxVtxRadius[0]);
3acc14d5 259 printf(" fMinDist2Vtxs: %f\n",fMinDist2Vtxs);
3acc14d5 260 printf("=======================================================\n");
261}
262//-------------------------------------------------------------------------