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