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