1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 "AliHLTITSVertexerZ.h"
21 #include "AliITSLoader.h"
22 #include <AliITSgeom.h>
23 #include <AliITSRecPoint.h>
24 #include <AliITSclusterV2.h>
26 //-------------------------------------------------------------------------
27 // Implementation of the HLT ITS vertexer class
29 // Origin: Cvetan Cheshkov, CERN, Cvetan.Cheshkov@cern.ch
30 //-------------------------------------------------------------------------
32 ClassImp(AliHLTITSVertexerZ)
34 AliHLTITSVertexerZ::AliHLTITSVertexerZ():AliITSVertexerZ(){
35 // Constructor in case that there is no runloader
40 SetFirstLayerModules();
41 SetSecondLayerModules();
55 AliHLTITSVertexerZ::AliHLTITSVertexerZ(TString filename,Float_t x0, Float_t y0):AliITSVertexerZ(filename,x0,y0)
57 // Standard Constructor
60 //______________________________________________________________________
61 AliESDVertex* AliHLTITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
62 // Defines the AliESDVertex for the current event
65 AliRunLoader *rl =AliRunLoader::GetRunLoader();
66 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
67 itsLoader->LoadRecPoints();
68 rl->GetEvent(evnumber);
71 fITS = (AliITS*)gAlice->GetModule("ITS");
73 Error("FindVertexForCurrentEvent","AliITS object was not found");
74 return fCurrentVertex;
79 // fITS->SetTreeAddress();
82 // AliITSgeom *geom = fITS->GetITSgeom();
83 AliITSgeom *geom = (AliITSgeom*)gDirectory->Get("AliITSgeom");
84 TTree *tR = itsLoader->TreeR();
86 return FindVertexForCurrentEvent(geom,tR);
89 //______________________________________________________________________
90 AliESDVertex* AliHLTITSVertexerZ::FindVertexForCurrentEvent(AliITSgeom *geom,TTree *tR){
91 // Defines the AliESDVertex for the current event
95 Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
96 Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
97 Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
98 Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
100 TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
102 branch = tR->GetBranch("Clusters");
103 branch->SetAddress(&clusters);
107 for(Int_t module= fFirstL1; module<=fLastL1;module++){
108 if(module%4==0 || module%4==3)continue;
109 // cout<<"Procesing module "<<module<<" ";
110 tR->GetEvent(module);
111 // cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
112 nrpL1+= clusters->GetEntriesFast();
114 for(Int_t module= fFirstL2; module<=fLastL2;module++){
115 tR->GetEvent(module);
116 nrpL2+= clusters->GetEntriesFast();
118 if(nrpL1 == 0 || nrpL2 == 0){
120 return fCurrentVertex;
123 Int_t nPhiBins = (Int_t)(TMath::Pi()/fDiffPhiMax);
124 Float_t phiBinSize = 2*TMath::Pi()/(Float_t)nPhiBins;
126 Int_t maxind1 = 2*nrpL1/nPhiBins;
127 Float_t **zc1 = new Float_t *[nPhiBins];
128 Float_t **phi1 = new Float_t *[nPhiBins];
129 Float_t **r1 = new Float_t *[nPhiBins];
130 Int_t *ind1 = new Int_t [nPhiBins];
131 Int_t maxind2 = 2*nrpL2/nPhiBins;
132 Float_t **zc2 = new Float_t *[nPhiBins];
133 Float_t **phi2 = new Float_t *[nPhiBins];
134 Float_t **r2 = new Float_t *[nPhiBins];
135 Int_t *ind2 = new Int_t [nPhiBins];
136 for(Int_t i=0;i<nPhiBins;i++) {
137 zc1[i] = new Float_t [maxind1];
138 phi1[i] = new Float_t [maxind1];
139 r1[i] = new Float_t [maxind1];
140 zc2[i] = new Float_t [maxind2];
141 phi2[i] = new Float_t [maxind2];
142 r2[i] = new Float_t [maxind2];
146 Float_t zshift[4] = {-10.708000, -3.536000, 3.536000, 10.708000};
149 memset(ind1,0,nPhiBins*sizeof(Int_t));
150 for(Int_t module= fFirstL1; module<=fLastL1;module++){
151 if(module%4==0 || module%4==3)continue;
152 tR->GetEvent(module);
153 Int_t nrecp1 = clusters->GetEntriesFast();
154 for(Int_t j=0;j<nrecp1;j++){
155 AliITSclusterV2 *recp = (AliITSclusterV2*)clusters->UncheckedAt(j);
156 lc[0]=-recp->GetY()+yshift;
157 lc[2]=-recp->GetZ()+zshift[module%4];
158 geom->LtoG(module,lc,gc);
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);
167 if(bin>=nPhiBins || bin<0) bin = 0;
168 Int_t ind = ind1[bin];
170 phi1[bin][ind] = phi;
171 zc1[bin][ind]=gc[2]/fStepFine;
172 r1[bin][ind]=sqrt(xc1*xc1+yc1*yc1);
179 memset(ind2,0,nPhiBins*sizeof(Int_t));
180 for(Int_t module= fFirstL2; module<=fLastL2;module++){
181 tR->GetEvent(module);
182 Int_t nrecp2 = clusters->GetEntriesFast();
183 for(Int_t j=0;j<nrecp2;j++){
184 AliITSclusterV2 *recp = (AliITSclusterV2*)clusters->UncheckedAt(j);
185 lc[0]=recp->GetY()+yshift;
186 lc[2]=-recp->GetZ()+zshift[module%4];
187 geom->LtoG(module,lc,gc);
193 Float_t phi = TMath::ATan2(gc[1],gc[0]);
194 if(phi<0)phi=2*TMath::Pi()+phi;
195 Int_t bin = (Int_t)(phi/phiBinSize+0.5);
196 if(bin>=nPhiBins || bin<0) bin = 0;
197 Int_t ind = ind2[bin];
199 phi2[bin][ind] = phi;
200 zc2[bin][ind]=gc[2]/fStepFine;
201 r2[bin][ind]=sqrt(xc2*xc2+yc2*yc2);
207 Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
208 Float_t lowz = fLowLim/fStepFine;
209 Int_t *harray = new Int_t[nbinfine];
210 memset(harray,0,nbinfine*sizeof(Int_t));
211 for(Int_t ibin=0;ibin<nPhiBins;ibin++) {
212 Float_t *pphi1 = phi1[ibin];
213 Float_t *pr1 = r1[ibin];
214 Float_t *pzc1 = zc1[ibin];
215 for(Int_t i=0;i<ind1[ibin];i++){
216 for(Int_t jbin=ibin;jbin<=(ibin+1);jbin++) {
218 if(jbin==nPhiBins) jbin2 = 0;
219 Float_t *pphi2 = phi2[jbin2];
220 Float_t *pr2 = r2[jbin2];
221 Float_t *pzc2 = zc2[jbin2];
222 for(Int_t j=0;j<ind2[jbin2];j++){
223 Float_t diff = TMath::Abs(pphi2[j]-pphi1[i]);
224 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
225 if(diff<fDiffPhiMax){
226 Float_t zr0=(pr2[j]*pzc1[i]-pr1[i]*pzc2[j])/(pr2[j]-pr1[i]);
227 Int_t bin = (Int_t)(zr0-lowz);
228 if(bin>=0 && bin<nbinfine){
236 for(Int_t i=0;i<nPhiBins;i++) {
253 Int_t nbinratio = (Int_t)(fStepCoarse/fStepFine+0.5);
254 Int_t nbincoarse = nbinfine/nbinratio;
256 if(fZCombc)delete fZCombc;
257 fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
258 if(fZCombf)delete fZCombf;
259 fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
262 for(Int_t bin=0;bin<nbinfine;bin++) {
263 fZCombf->SetBinContent(bin+1,(Stat_t)harray[bin]);
264 fZCombf->SetBinError(bin+1,TMath::Sqrt((Stat_t)harray[bin]));
265 contents+=(Stat_t)harray[bin];
267 if(counter==nbinratio) {
268 Int_t binc=bin/nbinratio;
269 fZCombc->SetBinContent(binc+1,contents);
270 fZCombc->SetBinError(binc+1,TMath::Sqrt(contents));
278 if(fZCombc->GetEntries()==0){
279 Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
281 return fCurrentVertex;
284 // cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
286 Int_t bi = fZCombc->GetMaximumBin();
287 Float_t centre = fZCombc->GetBinCenter(bi);
288 Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
289 Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
297 for(Int_t n=n1;n<=n2;n++){
298 fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
299 num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
300 fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
306 Float_t radi = fZsig/(num-1)-fZFound*fZFound/num/(num-1);
307 if(radi>0.)fZsig=TMath::Sqrt(radi);
311 goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
312 n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
313 n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
317 Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
320 // cout<<"Numer of Iterations "<<niter<<endl<<endl;
321 fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
322 fCurrentVertex->SetTitle("vertexer: HLT");
325 return fCurrentVertex;