adding map include file to solve compilation issue (bug https://savannah.cern.ch...
[u/mrichter/AliRoot.git] / HLT / ITS / AliHLTITSVertexerZ.cxx
CommitLineData
d6e2c707 1// $Id$
2
dd36288a 3/**************************************************************************
4 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * *
6 * Author: The ALICE Off-line Project. *
7 * Contributors are mentioned in the code where appropriate. *
8 * *
9 * Permission to use, copy, modify and distribute this software and its *
10 * documentation strictly for non-commercial purposes is hereby granted *
11 * without fee, provided that the above copyright notice appears in all *
12 * copies and that both the copyright notice and this permission notice *
13 * appear in the supporting documentation. The authors make no claims *
14 * about the suitability of this software for any purpose. It is *
15 * provided "as is" without express or implied warranty. *
16 **************************************************************************/
4aa41877 17#include "AliHLTITSVertexerZ.h"
dd36288a 18#include<TString.h>
19#include<TH1.h>
20#include<TMath.h>
21#include <AliRun.h>
22#include <AliITS.h>
23#include "AliITSLoader.h"
24#include <AliITSgeom.h>
25#include <AliITSRecPoint.h>
3136129d 26#include <AliITSclusterV2.h>
dd36288a 27
28//-------------------------------------------------------------------------
29// Implementation of the HLT ITS vertexer class
30//
31// Origin: Cvetan Cheshkov, CERN, Cvetan.Cheshkov@cern.ch
32//-------------------------------------------------------------------------
33
4aa41877 34ClassImp(AliHLTITSVertexerZ)
dd36288a 35
8b40a6b8 36AliHLTITSVertexerZ::AliHLTITSVertexerZ():
37 AliITSVertexerZ(),
38 fZCombf(0),
39 fStepFine(0)
40{
3136129d 41 // Constructor in case that there is no runloader
3136129d 42 SetBinWidthFine();
3136129d 43}
44
308c2f7c 45AliHLTITSVertexerZ::AliHLTITSVertexerZ(Float_t x0, Float_t y0):
46 AliITSVertexerZ(x0,y0),
8b40a6b8 47 fZCombf(0),
48 fStepFine(0)
dd36288a 49{
50 // Standard Constructor
8b40a6b8 51 SetBinWidthFine();
52}
53
54AliHLTITSVertexerZ::~AliHLTITSVertexerZ()
55{
56 // Destructor
57 if (fZCombf) delete fZCombf;
dd36288a 58}
59
60//______________________________________________________________________
4aa41877 61AliESDVertex* AliHLTITSVertexerZ::FindVertexForCurrentEvent(AliITSgeom *geom,TTree *tR){
3136129d 62 // Defines the AliESDVertex for the current event
63
64 fCurrentVertex = 0;
65
dd36288a 66 Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
67 Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
68 Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
69 Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
70
dd36288a 71 TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
72 TBranch *branch;
3136129d 73 branch = tR->GetBranch("Clusters");
74 branch->SetAddress(&clusters);
dd36288a 75
76 Int_t nrpL1 = 0;
77 Int_t nrpL2 = 0;
78 for(Int_t module= fFirstL1; module<=fLastL1;module++){
79 if(module%4==0 || module%4==3)continue;
80 // cout<<"Procesing module "<<module<<" ";
3136129d 81 tR->GetEvent(module);
dd36288a 82 // cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
3136129d 83 nrpL1+= clusters->GetEntriesFast();
dd36288a 84 }
85 for(Int_t module= fFirstL2; module<=fLastL2;module++){
3136129d 86 tR->GetEvent(module);
87 nrpL2+= clusters->GetEntriesFast();
dd36288a 88 }
89 if(nrpL1 == 0 || nrpL2 == 0){
90 ResetHistograms();
91 return fCurrentVertex;
92 }
93
94 Int_t nPhiBins = (Int_t)(TMath::Pi()/fDiffPhiMax);
95 Float_t phiBinSize = 2*TMath::Pi()/(Float_t)nPhiBins;
96
97 Int_t maxind1 = 2*nrpL1/nPhiBins;
98 Float_t **zc1 = new Float_t *[nPhiBins];
99 Float_t **phi1 = new Float_t *[nPhiBins];
100 Float_t **r1 = new Float_t *[nPhiBins];
101 Int_t *ind1 = new Int_t [nPhiBins];
102 Int_t maxind2 = 2*nrpL2/nPhiBins;
103 Float_t **zc2 = new Float_t *[nPhiBins];
104 Float_t **phi2 = new Float_t *[nPhiBins];
105 Float_t **r2 = new Float_t *[nPhiBins];
106 Int_t *ind2 = new Int_t [nPhiBins];
107 for(Int_t i=0;i<nPhiBins;i++) {
108 zc1[i] = new Float_t [maxind1];
109 phi1[i] = new Float_t [maxind1];
110 r1[i] = new Float_t [maxind1];
111 zc2[i] = new Float_t [maxind2];
112 phi2[i] = new Float_t [maxind2];
113 r2[i] = new Float_t [maxind2];
114 }
3136129d 115
116 Float_t yshift = 0;
117 Float_t zshift[4] = {-10.708000, -3.536000, 3.536000, 10.708000};
dd36288a 118
3136129d 119 yshift = 0.248499;
dd36288a 120 memset(ind1,0,nPhiBins*sizeof(Int_t));
121 for(Int_t module= fFirstL1; module<=fLastL1;module++){
122 if(module%4==0 || module%4==3)continue;
3136129d 123 tR->GetEvent(module);
124 Int_t nrecp1 = clusters->GetEntriesFast();
dd36288a 125 for(Int_t j=0;j<nrecp1;j++){
3136129d 126 AliITSclusterV2 *recp = (AliITSclusterV2*)clusters->UncheckedAt(j);
127 lc[0]=-recp->GetY()+yshift;
128 lc[2]=-recp->GetZ()+zshift[module%4];
dd36288a 129 geom->LtoG(module,lc,gc);
308c2f7c 130 gc[0]-=GetNominalPos()[0];
131 gc[1]-=GetNominalPos()[1];
dd36288a 132 Float_t xc1,yc1;
133 xc1=gc[0];
134 yc1=gc[1];
135 Float_t phi = TMath::ATan2(gc[1],gc[0]);
136 if(phi<0)phi=2*TMath::Pi()+phi;
137 Int_t bin = (Int_t)(phi/phiBinSize);
3136129d 138 if(bin>=nPhiBins || bin<0) bin = 0;
dd36288a 139 Int_t ind = ind1[bin];
140 if(ind<maxind1) {
141 phi1[bin][ind] = phi;
142 zc1[bin][ind]=gc[2]/fStepFine;
143 r1[bin][ind]=sqrt(xc1*xc1+yc1*yc1);
144 ind1[bin]++;
145 }
146 }
3136129d 147 clusters->Delete();
dd36288a 148 }
3136129d 149 yshift = 3.096207;
dd36288a 150 memset(ind2,0,nPhiBins*sizeof(Int_t));
151 for(Int_t module= fFirstL2; module<=fLastL2;module++){
3136129d 152 tR->GetEvent(module);
153 Int_t nrecp2 = clusters->GetEntriesFast();
dd36288a 154 for(Int_t j=0;j<nrecp2;j++){
3136129d 155 AliITSclusterV2 *recp = (AliITSclusterV2*)clusters->UncheckedAt(j);
156 lc[0]=recp->GetY()+yshift;
157 lc[2]=-recp->GetZ()+zshift[module%4];
dd36288a 158 geom->LtoG(module,lc,gc);
308c2f7c 159 gc[0]-=GetNominalPos()[0];
160 gc[1]-=GetNominalPos()[1];
dd36288a 161 Float_t xc2,yc2;
162 xc2=gc[0];
163 yc2=gc[1];
164 Float_t phi = TMath::ATan2(gc[1],gc[0]);
165 if(phi<0)phi=2*TMath::Pi()+phi;
166 Int_t bin = (Int_t)(phi/phiBinSize+0.5);
3136129d 167 if(bin>=nPhiBins || bin<0) bin = 0;
dd36288a 168 Int_t ind = ind2[bin];
169 if(ind<maxind2) {
170 phi2[bin][ind] = phi;
171 zc2[bin][ind]=gc[2]/fStepFine;
172 r2[bin][ind]=sqrt(xc2*xc2+yc2*yc2);
173 ind2[bin]++;
174 }
175 }
3136129d 176 clusters->Delete();
dd36288a 177 }
178 Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
179 Float_t lowz = fLowLim/fStepFine;
dd36288a 180 Int_t *harray = new Int_t[nbinfine];
181 memset(harray,0,nbinfine*sizeof(Int_t));
182 for(Int_t ibin=0;ibin<nPhiBins;ibin++) {
183 Float_t *pphi1 = phi1[ibin];
184 Float_t *pr1 = r1[ibin];
185 Float_t *pzc1 = zc1[ibin];
186 for(Int_t i=0;i<ind1[ibin];i++){
187 for(Int_t jbin=ibin;jbin<=(ibin+1);jbin++) {
188 Int_t jbin2 = jbin;
189 if(jbin==nPhiBins) jbin2 = 0;
190 Float_t *pphi2 = phi2[jbin2];
191 Float_t *pr2 = r2[jbin2];
192 Float_t *pzc2 = zc2[jbin2];
193 for(Int_t j=0;j<ind2[jbin2];j++){
194 Float_t diff = TMath::Abs(pphi2[j]-pphi1[i]);
195 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
196 if(diff<fDiffPhiMax){
197 Float_t zr0=(pr2[j]*pzc1[i]-pr1[i]*pzc2[j])/(pr2[j]-pr1[i]);
3136129d 198 Int_t bin = (Int_t)(zr0-lowz);
199 if(bin>=0 && bin<nbinfine){
dd36288a 200 harray[bin]++;
201 }
202 }
203 }
204 }
205 }
206 }
207 for(Int_t i=0;i<nPhiBins;i++) {
208 delete [] zc1[i];
209 delete [] phi1[i];
210 delete [] r1[i];
211 delete [] zc2[i];
212 delete [] phi2[i];
213 delete [] r2[i];
214 }
215 delete [] zc1;
216 delete [] phi1;
217 delete [] r1;
218 delete [] ind1;
219 delete [] zc2;
220 delete [] phi2;
221 delete [] r2;
222 delete [] ind2;
223
224 Int_t nbinratio = (Int_t)(fStepCoarse/fStepFine+0.5);
225 Int_t nbincoarse = nbinfine/nbinratio;
226
227 if(fZCombc)delete fZCombc;
228 fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
229 if(fZCombf)delete fZCombf;
230 fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
231 Stat_t contents=0;
232 Int_t counter=0;
233 for(Int_t bin=0;bin<nbinfine;bin++) {
234 fZCombf->SetBinContent(bin+1,(Stat_t)harray[bin]);
235 fZCombf->SetBinError(bin+1,TMath::Sqrt((Stat_t)harray[bin]));
236 contents+=(Stat_t)harray[bin];
237 counter++;
238 if(counter==nbinratio) {
239 Int_t binc=bin/nbinratio;
240 fZCombc->SetBinContent(binc+1,contents);
241 fZCombc->SetBinError(binc+1,TMath::Sqrt(contents));
242 contents=0;
243 counter=0;
244 }
245 }
246
247 delete [] harray;
248
249 if(fZCombc->GetEntries()==0){
250 Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
251 ResetHistograms();
252 return fCurrentVertex;
253 }
254 // else {
255 // cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
256 // }
257 Int_t bi = fZCombc->GetMaximumBin();
258 Float_t centre = fZCombc->GetBinCenter(bi);
259 Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
260 Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
261 Int_t niter = 0;
262 Bool_t goon = kTRUE;
263 Int_t num;
264 while(goon){
265 fZFound = 0.;
266 fZsig = 0.;
267 num=0;
268 for(Int_t n=n1;n<=n2;n++){
269 fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
270 num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
271 fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
272 }
273 if(num<2){
274 fZsig = 0.;
275 }
276 else {
277 Float_t radi = fZsig/(num-1)-fZFound*fZFound/num/(num-1);
278 if(radi>0.)fZsig=TMath::Sqrt(radi);
279 else fZsig=0.;
280 fZFound/=num;
281 }
282 goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
283 n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
284 n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
285 niter++;
286 if(niter>=10){
287 goon = kFALSE;
288 Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
289 }
290 }
291 // cout<<"Numer of Iterations "<<niter<<endl<<endl;
292 fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
293 fCurrentVertex->SetTitle("vertexer: HLT");
294 ResetHistograms();
295 PrintStatus();
296 return fCurrentVertex;
297}