Added variable init to copy constructor
[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
28// vertex for cosmic events in which the muon crosses the SPD inner
29// layer (SPD1). 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 SPD1 and define the fake vertex as
33// the mid-point of the straight line joining the two clusters.
34// We reject the backgroung by requiring at least one clusters on SPD2
35// closer than fMaxDistOnSPD2 to the tracklet prolongation.
36// We can reject (potentially pathological) events with the muon track
37// tangential to the SPD1 layer by the requiring the radial position of
38// the vertex to be smaller than fMaxVtxRadius.
39// Due to background clusters, more than one vertex per event can
40// be found. We consider the first found.
41// The number of contributors set in the AliESDVertex object is the
42// number of vertices found in the event; if this number is <1,
43// the procedure could not find a vertex position and by default
44// the vertex coordinates are set to (0,0,0) with large errors (100,100,100)
45// Number of contributors = 0 --> No SPD1 tracklets matching criteria
46// Number of contributors = -1 --> No SPD1 recpoints
47//
48// Origin: A.Dainese, andrea.dainese@lnl.infn.it
49//-------------------------------------------------------------------------
50
51ClassImp(AliITSVertexerCosmics)
52
53//-------------------------------------------------------------------------
54AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
55fFirstSPD1(0),
56fLastSPD1(0),
57fFirstSPD2(0),
58fLastSPD2(0),
59fMaxDistOnSPD2(0),
60fMaxVtxRadius(0),
61fMinDist2Vtxs(0)
62{
63 // Default constructor
64 SetSPD1Modules();
65 SetSPD2Modules();
66 SetMaxDistOnSPD2();
67 SetMaxVtxRadius();
68 SetMinDist2Vtxs();
69}
70//--------------------------------------------------------------------------
71AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(Int_t evnumber)
72{
73 // Defines the AliESDVertex for the current event
74
75 fCurrentVertex = 0;
76 AliRunLoader *rl =AliRunLoader::GetRunLoader();
77 AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
3acc14d5 78 itsLoader->LoadRecPoints();
79 rl->GetEvent(evnumber);
80
81 TTree *rpTree = itsLoader->TreeR();
82
83 TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
84 rpTree->SetBranchAddress("ITSRecPoints",&recpoints);
85
86 Double_t xclspd1[100],yclspd1[100],zclspd1[100],modclspd1[100];
87 Int_t nclspd1stored=0;
88 Double_t xclspd2[100],yclspd2[100],zclspd2[100],modclspd2[100];
89 Int_t nclspd2stored=0;
90 Int_t nrecpoints,nrecpointsSPD1=0;
91
92 Double_t gc[3]={0.,0.,0.};
93 Double_t lc[3]={0.,0.,0.};
94 Int_t lay,lad,det;
95
96 Double_t x[100],y[100],z[100],p1[3],p2[3],p3[3];
97 Int_t nvtxs;
98 Bool_t good,matchtospd2;
99 Double_t xvtx,yvtx,zvtx,rvtx;
100
101 for(Int_t imodule=fFirstSPD1; imodule<fLastSPD2; imodule++) { // SPD
102 rpTree->GetEvent(imodule);
0333e459 103 AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
3acc14d5 104 nrecpoints=recpoints->GetEntriesFast();
105 if(imodule<fLastSPD1) nrecpointsSPD1 += nrecpoints;
106 //printf("cosmics: module %d clusters %d\n",imodule,nrecpoints);
107 for(Int_t irp=0; irp<nrecpoints; irp++) {
108 AliITSRecPoint *rp=(AliITSRecPoint*)recpoints->UncheckedAt(irp);
109 // Local coordinates of this recpoint
110 lc[0]=rp->GetDetLocalX();
111 lc[2]=rp->GetDetLocalZ();
0333e459 112 AliITSgeomTGeo::LocalToGlobal(imodule,lc,gc); // global coordinates
3acc14d5 113 if(lay==1) { // store SPD1 clusters
114 xclspd1[nclspd1stored]=gc[0];
115 yclspd1[nclspd1stored]=gc[1];
116 zclspd1[nclspd1stored]=gc[2];
117 modclspd1[nclspd1stored]=imodule;
118 nclspd1stored++;
119 }
120 if(lay==2) { // store SPD2 clusters
121 xclspd2[nclspd2stored]=gc[0];
122 yclspd2[nclspd2stored]=gc[1];
123 zclspd2[nclspd2stored]=gc[2];
124 modclspd2[nclspd2stored]=imodule;
125 nclspd2stored++;
126 }
127 if(nclspd1stored>100 || nclspd2stored>100)
128 AliFatal("More than 100 clusters per layer in SPD");
129 }// end clusters in a module
130 }// end SPD modules for a given event
131
132 // build fake vertices
133 nvtxs=0;
134 // SPD1 - first cluster
135 for(Int_t i1spd1=0; i1spd1<nclspd1stored; i1spd1++) {
136 p1[0]=xclspd1[i1spd1]; p1[1]=yclspd1[i1spd1]; p1[2]=zclspd1[i1spd1];
137 // SPD1 - second cluster
138 for(Int_t i2spd1=i1spd1+1; i2spd1<nclspd1stored; i2spd1++) {
139 if(modclspd1[i1spd1]==modclspd1[i2spd1]) continue;
140 p2[0]=xclspd1[i2spd1]; p2[1]=yclspd1[i2spd1]; p2[2]=zclspd1[i2spd1];
141 // look for point on SPD2
142 AliStrLine spd1line(p1,p2,kTRUE);
143 matchtospd2 = kFALSE;
144 for(Int_t ispd2=0; ispd2<nclspd2stored; ispd2++) {
145 p3[0]=xclspd2[ispd2]; p3[1]=yclspd2[ispd2]; p3[2]=zclspd2[ispd2];
146 //printf(" %f\n",spd1line.GetDistFromPoint(p3));
147 if(spd1line.GetDistFromPoint(p3)<fMaxDistOnSPD2)
148 { matchtospd2 = kTRUE; break; }
149 }
150 if(!matchtospd2) continue;
151 xvtx = 0.5*(xclspd1[i1spd1]+xclspd1[i2spd1]);
152 yvtx = 0.5*(yclspd1[i1spd1]+yclspd1[i2spd1]);
153 zvtx = 0.5*(zclspd1[i1spd1]+zclspd1[i2spd1]);
154 rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
155 if(rvtx>fMaxVtxRadius) continue;
156 good = kTRUE;
157 for(Int_t iv=0; iv<nvtxs; iv++) {
158 if(TMath::Sqrt((xvtx- x[iv])*(xvtx- x[iv])+
159 (yvtx- y[iv])*(yvtx- y[iv])+
160 (zvtx- z[iv])*(zvtx- z[iv])) < fMinDist2Vtxs)
161 good = kFALSE;
162 }
163 if(good) {
164 x[nvtxs]=xvtx;
165 y[nvtxs]=yvtx;
166 z[nvtxs]=zvtx;
167 nvtxs++;
168 }
169 } // SPD1 - second cluster
170 } // SPD1 - first cluster
171
172
173 Double_t pos[3]={0.,0.,0.};
174 Double_t err[3]={100.,100.,100.};
175 if(nvtxs) {
176 pos[0]=x[0];
177 pos[1]=y[0];
178 pos[2]=z[0];
179 err[0]=0.1;
180 err[1]=0.1;
181 err[2]=0.1;
182 }
183 if(!nrecpointsSPD1) nvtxs=-1;
184 fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
185 fCurrentVertex->SetNContributors(nvtxs);
186 fCurrentVertex->SetTitle("cosmics fake vertex");
187
188 //if(nvtxs>0) fCurrentVertex->Print();
189
190 delete recpoints;
191 itsLoader->UnloadRecPoints();
192
193 return fCurrentVertex;
194}
195//-------------------------------------------------------------------------
196void AliITSVertexerCosmics::FindVertices()
197{
198 // computes the vertices of the events in the range FirstEvent - LastEvent
199 AliRunLoader *rl = AliRunLoader::GetRunLoader();
200 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
201 itsLoader->ReloadRecPoints();
202 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
203 // printf("Processing event %d\n",i);
204 rl->GetEvent(i);
205 FindVertexForCurrentEvent(i);
206 if(fCurrentVertex){
207 WriteCurrentVertex();
208 }
209 }
210}
211//-------------------------------------------------------------------------
212void AliITSVertexerCosmics::PrintStatus() const
213{
214 // Print current status
215 printf("=======================================================\n");
216 printf(" fMaxDistOnSPD2: %f\n",fMaxDistOnSPD2);
217 printf(" fMaxVtxRadius: %f\n",fMaxVtxRadius);
218 printf(" fMinDist2Vtxs: %f\n",fMinDist2Vtxs);
219 printf(" First layer first and last modules: %d, %d\n",fFirstSPD1,fLastSPD1);
220 printf(" Second layer first and last modules: %d, %d\n",fFirstSPD2,fLastSPD2);
221 printf("=======================================================\n");
222}
223//-------------------------------------------------------------------------