merging RecPoints and ClustersV2. All ClusterFinders produce AliITSRecPoints objects...
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerZ.cxx
CommitLineData
0c6af5c9 1/**************************************************************************
2 * Copyright(c) 1998-2003, 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#include <AliITSVertexerZ.h>
0c6af5c9 16#include <TString.h>
0c6af5c9 17#include "AliITSLoader.h"
0c6af5c9 18#include<TBranch.h>
19#include<TClonesArray.h>
20#include<TH1.h>
0c6af5c9 21#include<TTree.h>
22#include <AliITSgeom.h>
7d62fb64 23#include "AliITSDetTypeRec.h"
0c6af5c9 24#include <AliITSRecPoint.h>
25
26/////////////////////////////////////////////////////////////////
27// this class implements a fast method to determine
28// the Z coordinate of the primary vertex
29// for p-p collisions it seems to give comparable or better results
30// with respect to what obtained with AliITSVertexerPPZ
31// It can be used successfully with Pb-Pb collisions
32////////////////////////////////////////////////////////////////
33
34ClassImp(AliITSVertexerZ)
35
36
37
38//______________________________________________________________________
39AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer() {
40 // Default constructor
41 SetDiffPhiMax(0);
42 fX0 = 0.;
43 fY0 = 0.;
44 SetFirstLayerModules(0);
45 SetSecondLayerModules(0);
46 fZFound = 0;
47 fZsig = 0.;
0c6af5c9 48 fZCombc = 0;
49 fZCombf = 0;
50 SetLowLimit(0.);
51 SetHighLimit(0.);
52 SetBinWidthCoarse(0.);
53 SetBinWidthFine(0.);
54 SetTolerance(0.);
55}
56
57//______________________________________________________________________
58AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
59 // Standard Constructor
60 SetDiffPhiMax();
61 fX0 = x0;
62 fY0 = y0;
63 SetFirstLayerModules();
64 SetSecondLayerModules();
65 fZFound = 0;
66 fZsig = 0.;
0c6af5c9 67 fZCombc = 0;
68 fZCombf = 0;
69 SetLowLimit();
70 SetHighLimit();
71 SetBinWidthCoarse();
72 SetBinWidthFine();
73 SetTolerance();
74
75}
76
77//______________________________________________________________________
78AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr) {
79 // Copy constructor
80 // Copies are not allowed. The method is protected to avoid misuse.
81 Error("AliITSVertexerZ","Copy constructor not allowed\n");
82}
83
84//______________________________________________________________________
85AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ& /* vtxr */){
86 // Assignment operator
87 // Assignment is not allowed. The method is protected to avoid misuse.
88 Error("= operator","Assignment operator not allowed\n");
89 return *this;
90}
91
92
93//______________________________________________________________________
94AliITSVertexerZ::~AliITSVertexerZ() {
95 // Default Destructor
7d62fb64 96 //fITS = 0;
0c6af5c9 97 if(fZCombc)delete fZCombc;
98 if(fZCombf)delete fZCombf;
99}
100
101//______________________________________________________________________
d681bb2d 102AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
103 // Defines the AliESDVertex for the current event
7d62fb64 104
0c6af5c9 105 fCurrentVertex = 0;
106 AliRunLoader *rl =AliRunLoader::GetRunLoader();
7d62fb64 107 AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
108 TDirectory * olddir = gDirectory;
109 rl->CdGAFile();
110 AliITSgeom* geom = (AliITSgeom*)gDirectory->Get("AliITSgeom");
111 olddir->cd();
112
2257f27e 113 itsLoader->LoadRecPoints();
0c6af5c9 114 rl->GetEvent(evnumber);
115
7d62fb64 116 AliITSDetTypeRec detTypeRec;
0c6af5c9 117
118 TTree *tR = itsLoader->TreeR();
7d62fb64 119 detTypeRec.SetTreeAddressR(tR);
0c6af5c9 120 TClonesArray *itsRec = 0;
121 Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
122 Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
123 Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
124 Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
125
7d62fb64 126 itsRec = detTypeRec.RecPoints();
00a7cc50 127 // TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
2257f27e 128 TBranch *branch;
00a7cc50 129 //if(fUseV2Clusters){
130 // branch = tR->GetBranch("ITSRecPoints");
131 // branch->SetAddress(&clusters);
132 //}
133 //else {
134 branch = tR->GetBranch("ITSRecPoints");
135 //}
0c6af5c9 136
137 Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
138 Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
139 if(fZCombc)delete fZCombc;
140 fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
141 if(fZCombf)delete fZCombf;
142 fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
143
144 Int_t nrpL1 = 0;
145 Int_t nrpL2 = 0;
146 for(Int_t module= fFirstL1; module<=fLastL1;module++){
147 if(module%4==0 || module%4==3)continue;
2257f27e 148 // cout<<"Procesing module "<<module<<" ";
0c6af5c9 149 branch->GetEvent(module);
2257f27e 150 // cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
00a7cc50 151 //if(fUseV2Clusters){
152 // Clusters2RecPoints(clusters,module,itsRec);
153 // }
0c6af5c9 154 nrpL1+= itsRec->GetEntries();
7d62fb64 155 detTypeRec.ResetRecPoints();
0c6af5c9 156 }
157 for(Int_t module= fFirstL2; module<=fLastL2;module++){
158 branch->GetEvent(module);
00a7cc50 159 //if(fUseV2Clusters){
160 // Clusters2RecPoints(clusters,module,itsRec);
161 //}
0c6af5c9 162 nrpL2+= itsRec->GetEntries();
7d62fb64 163 detTypeRec.ResetRecPoints();
0c6af5c9 164 }
0c6af5c9 165 if(nrpL1 == 0 || nrpL2 == 0){
0c6af5c9 166 ResetHistograms();
167 return fCurrentVertex;
168 }
169 Float_t *xc1 = new Float_t [nrpL1];
170 Float_t *yc1 = new Float_t [nrpL1];
171 Float_t *zc1 = new Float_t [nrpL1];
172 Float_t *phi1 = new Float_t [nrpL1];
173 Float_t *xc2 = new Float_t [nrpL2];
174 Float_t *yc2 = new Float_t [nrpL2];
175 Float_t *zc2 = new Float_t [nrpL2];
176 Float_t *phi2 = new Float_t [nrpL2];
177 Int_t ind = 0;
178 for(Int_t module= fFirstL1; module<=fLastL1;module++){
179 if(module%4==0 || module%4==3)continue;
180 branch->GetEvent(module);
00a7cc50 181 //if(fUseV2Clusters){
182 // Clusters2RecPoints(clusters,module,itsRec);
183 //}
0c6af5c9 184 Int_t nrecp1 = itsRec->GetEntries();
185 for(Int_t j=0;j<nrecp1;j++){
186 AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
00a7cc50 187 lc[0]=recp->GetDetLocalX();
188 lc[2]=recp->GetDetLocalZ();
0c6af5c9 189 geom->LtoG(module,lc,gc);
190 gc[0]-=fX0;
191 gc[1]-=fY0;
192 xc1[ind]=gc[0];
193 yc1[ind]=gc[1];
194 zc1[ind]=gc[2];
195 phi1[ind] = TMath::ATan2(gc[1],gc[0]);
196 if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
197 ind++;
198 }
7d62fb64 199 detTypeRec.ResetRecPoints();
0c6af5c9 200 }
201 ind = 0;
202 for(Int_t module= fFirstL2; module<=fLastL2;module++){
203 branch->GetEvent(module);
00a7cc50 204 //if(fUseV2Clusters){
205 // Clusters2RecPoints(clusters,module,itsRec);
206 //}
0c6af5c9 207 Int_t nrecp2 = itsRec->GetEntries();
208 for(Int_t j=0;j<nrecp2;j++){
209 AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
00a7cc50 210 lc[0]=recp->GetDetLocalX();
211 lc[2]=recp->GetDetLocalZ();
0c6af5c9 212 geom->LtoG(module,lc,gc);
213 gc[0]-=fX0;
214 gc[1]-=fY0;
215 xc2[ind]=gc[0];
216 yc2[ind]=gc[1];
217 zc2[ind]=gc[2];
218 phi2[ind] = TMath::ATan2(gc[1],gc[0]);
219 if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
220 ind++;
221 }
7d62fb64 222 detTypeRec.ResetRecPoints();
0c6af5c9 223 }
224 for(Int_t i=0;i<nrpL1;i++){
225 Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]);
226 for(Int_t j=0;j<nrpL2;j++){
227 Float_t diff = TMath::Abs(phi2[j]-phi1[i]);
228 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
229 if(diff<fDiffPhiMax){
230 Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]);
231 Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1);
232 fZCombf->Fill(zr0);
233 fZCombc->Fill(zr0);
234 }
235 }
236 }
237 delete [] xc1;
238 delete [] yc1;
239 delete [] zc1;
240 delete [] phi1;
241 delete [] xc2;
242 delete [] yc2;
243 delete [] zc2;
244 delete [] phi2;
245 if(fZCombc->GetEntries()==0){
246 Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
247 ResetHistograms();
248 return fCurrentVertex;
249 }
2257f27e 250 // else {
251 // cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
252 // }
0c6af5c9 253 Int_t bi = fZCombc->GetMaximumBin();
254 Float_t centre = fZCombc->GetBinCenter(bi);
255 Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
256 Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
257 Int_t niter = 0;
258 Bool_t goon = kTRUE;
259 Int_t num;
260 while(goon){
261 fZFound = 0.;
262 fZsig = 0.;
263 num=0;
264 for(Int_t n=n1;n<=n2;n++){
265 fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
266 num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
267 fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
268 }
269 if(num<2){
0c6af5c9 270 fZsig = 0.;
271 }
272 else {
979e3647 273 Float_t radi = fZsig/(num-1)-fZFound*fZFound/num/(num-1);
274 if(radi>0.)fZsig=TMath::Sqrt(radi);
275 else fZsig=0.;
0c6af5c9 276 fZFound/=num;
277 }
278 goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
0c6af5c9 279 n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
280 n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
281 niter++;
282 if(niter>=10){
283 goon = kFALSE;
284 Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
285 }
286 }
2257f27e 287 // cout<<"Numer of Iterations "<<niter<<endl<<endl;
d681bb2d 288 fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
0c6af5c9 289 fCurrentVertex->SetTitle("vertexer: B");
290 ResetHistograms();
291 return fCurrentVertex;
292}
293
294//_____________________________________________________________________
295void AliITSVertexerZ::ResetHistograms(){
296 // delete TH1 data members
297 if(fZCombc)delete fZCombc;
298 if(fZCombf)delete fZCombf;
299 fZCombc = 0;
300 fZCombf = 0;
301}
302
303//______________________________________________________________________
304void AliITSVertexerZ::FindVertices(){
305 // computes the vertices of the events in the range FirstEvent - LastEvent
306 AliRunLoader *rl = AliRunLoader::GetRunLoader();
307 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
308 itsLoader->ReloadRecPoints();
309 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
310 cout<<"Processing event "<<i<<endl;
311 rl->GetEvent(i);
312 FindVertexForCurrentEvent(i);
313 if(fCurrentVertex){
314 WriteCurrentVertex();
315 }
316 else {
317 if(fDebug>0){
318 cout<<"Vertex not found for event "<<i<<endl;
319 cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
320 }
321 }
322 }
323}
324
325//________________________________________________________
326void AliITSVertexerZ::PrintStatus() const {
327 // Print current status
328 cout <<"=======================================================\n";
329 cout <<" First layer first and last modules: "<<fFirstL1<<", ";
330 cout <<fLastL1<<endl;
331 cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
332 cout <<fLastL2<<endl;
333 cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
334 cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
335 cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
336 cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
337 cout <<" Debug flag: "<<fDebug<<endl;
338 cout <<"First event to be processed "<<fFirstEvent;
339 cout <<"\n Last event to be processed "<<fLastEvent<<endl;
340
341
342 cout <<"=======================================================\n";
343}
344