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