fix for accessing T0 data (C.Zampolli)
[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"
21#include "AliITSgeom.h"
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");
78 AliITSgeom* geom = itsLoader->GetITSgeom();
79 itsLoader->LoadRecPoints();
80 rl->GetEvent(evnumber);
81
82 TTree *rpTree = itsLoader->TreeR();
83
84 TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
85 rpTree->SetBranchAddress("ITSRecPoints",&recpoints);
86
87 Double_t xclspd1[100],yclspd1[100],zclspd1[100],modclspd1[100];
88 Int_t nclspd1stored=0;
89 Double_t xclspd2[100],yclspd2[100],zclspd2[100],modclspd2[100];
90 Int_t nclspd2stored=0;
91 Int_t nrecpoints,nrecpointsSPD1=0;
92
93 Double_t gc[3]={0.,0.,0.};
94 Double_t lc[3]={0.,0.,0.};
95 Int_t lay,lad,det;
96
97 Double_t x[100],y[100],z[100],p1[3],p2[3],p3[3];
98 Int_t nvtxs;
99 Bool_t good,matchtospd2;
100 Double_t xvtx,yvtx,zvtx,rvtx;
101
102 for(Int_t imodule=fFirstSPD1; imodule<fLastSPD2; imodule++) { // SPD
103 rpTree->GetEvent(imodule);
104 geom->GetModuleId(imodule,lay,lad,det);
105 nrecpoints=recpoints->GetEntriesFast();
106 if(imodule<fLastSPD1) nrecpointsSPD1 += nrecpoints;
107 //printf("cosmics: module %d clusters %d\n",imodule,nrecpoints);
108 for(Int_t irp=0; irp<nrecpoints; irp++) {
109 AliITSRecPoint *rp=(AliITSRecPoint*)recpoints->UncheckedAt(irp);
110 // Local coordinates of this recpoint
111 lc[0]=rp->GetDetLocalX();
112 lc[2]=rp->GetDetLocalZ();
113 geom->LtoG(imodule,lc,gc); // global coordinates
114 if(lay==1) { // store SPD1 clusters
115 xclspd1[nclspd1stored]=gc[0];
116 yclspd1[nclspd1stored]=gc[1];
117 zclspd1[nclspd1stored]=gc[2];
118 modclspd1[nclspd1stored]=imodule;
119 nclspd1stored++;
120 }
121 if(lay==2) { // store SPD2 clusters
122 xclspd2[nclspd2stored]=gc[0];
123 yclspd2[nclspd2stored]=gc[1];
124 zclspd2[nclspd2stored]=gc[2];
125 modclspd2[nclspd2stored]=imodule;
126 nclspd2stored++;
127 }
128 if(nclspd1stored>100 || nclspd2stored>100)
129 AliFatal("More than 100 clusters per layer in SPD");
130 }// end clusters in a module
131 }// end SPD modules for a given event
132
133 // build fake vertices
134 nvtxs=0;
135 // SPD1 - first cluster
136 for(Int_t i1spd1=0; i1spd1<nclspd1stored; i1spd1++) {
137 p1[0]=xclspd1[i1spd1]; p1[1]=yclspd1[i1spd1]; p1[2]=zclspd1[i1spd1];
138 // SPD1 - second cluster
139 for(Int_t i2spd1=i1spd1+1; i2spd1<nclspd1stored; i2spd1++) {
140 if(modclspd1[i1spd1]==modclspd1[i2spd1]) continue;
141 p2[0]=xclspd1[i2spd1]; p2[1]=yclspd1[i2spd1]; p2[2]=zclspd1[i2spd1];
142 // look for point on SPD2
143 AliStrLine spd1line(p1,p2,kTRUE);
144 matchtospd2 = kFALSE;
145 for(Int_t ispd2=0; ispd2<nclspd2stored; ispd2++) {
146 p3[0]=xclspd2[ispd2]; p3[1]=yclspd2[ispd2]; p3[2]=zclspd2[ispd2];
147 //printf(" %f\n",spd1line.GetDistFromPoint(p3));
148 if(spd1line.GetDistFromPoint(p3)<fMaxDistOnSPD2)
149 { matchtospd2 = kTRUE; break; }
150 }
151 if(!matchtospd2) continue;
152 xvtx = 0.5*(xclspd1[i1spd1]+xclspd1[i2spd1]);
153 yvtx = 0.5*(yclspd1[i1spd1]+yclspd1[i2spd1]);
154 zvtx = 0.5*(zclspd1[i1spd1]+zclspd1[i2spd1]);
155 rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
156 if(rvtx>fMaxVtxRadius) continue;
157 good = kTRUE;
158 for(Int_t iv=0; iv<nvtxs; iv++) {
159 if(TMath::Sqrt((xvtx- x[iv])*(xvtx- x[iv])+
160 (yvtx- y[iv])*(yvtx- y[iv])+
161 (zvtx- z[iv])*(zvtx- z[iv])) < fMinDist2Vtxs)
162 good = kFALSE;
163 }
164 if(good) {
165 x[nvtxs]=xvtx;
166 y[nvtxs]=yvtx;
167 z[nvtxs]=zvtx;
168 nvtxs++;
169 }
170 } // SPD1 - second cluster
171 } // SPD1 - first cluster
172
173
174 Double_t pos[3]={0.,0.,0.};
175 Double_t err[3]={100.,100.,100.};
176 if(nvtxs) {
177 pos[0]=x[0];
178 pos[1]=y[0];
179 pos[2]=z[0];
180 err[0]=0.1;
181 err[1]=0.1;
182 err[2]=0.1;
183 }
184 if(!nrecpointsSPD1) nvtxs=-1;
185 fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
186 fCurrentVertex->SetNContributors(nvtxs);
187 fCurrentVertex->SetTitle("cosmics fake vertex");
188
189 //if(nvtxs>0) fCurrentVertex->Print();
190
191 delete recpoints;
192 itsLoader->UnloadRecPoints();
193
194 return fCurrentVertex;
195}
196//-------------------------------------------------------------------------
197void AliITSVertexerCosmics::FindVertices()
198{
199 // computes the vertices of the events in the range FirstEvent - LastEvent
200 AliRunLoader *rl = AliRunLoader::GetRunLoader();
201 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
202 itsLoader->ReloadRecPoints();
203 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
204 // printf("Processing event %d\n",i);
205 rl->GetEvent(i);
206 FindVertexForCurrentEvent(i);
207 if(fCurrentVertex){
208 WriteCurrentVertex();
209 }
210 }
211}
212//-------------------------------------------------------------------------
213void AliITSVertexerCosmics::PrintStatus() const
214{
215 // Print current status
216 printf("=======================================================\n");
217 printf(" fMaxDistOnSPD2: %f\n",fMaxDistOnSPD2);
218 printf(" fMaxVtxRadius: %f\n",fMaxVtxRadius);
219 printf(" fMinDist2Vtxs: %f\n",fMinDist2Vtxs);
220 printf(" First layer first and last modules: %d, %d\n",fFirstSPD1,fLastSPD1);
221 printf(" Second layer first and last modules: %d, %d\n",fFirstSPD2,fLastSPD2);
222 printf("=======================================================\n");
223}
224//-------------------------------------------------------------------------