L3 becomes HLT
[u/mrichter/AliRoot.git] / HLT / ITS / AliHLTITSVertexerZ.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 "AliHLTITSVertexerZ.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 #include <AliITSclusterV2.h>
25
26 //-------------------------------------------------------------------------
27 //                Implementation of the HLT ITS vertexer class
28 //
29 //          Origin: Cvetan Cheshkov, CERN, Cvetan.Cheshkov@cern.ch
30 //-------------------------------------------------------------------------
31
32 ClassImp(AliHLTITSVertexerZ)
33
34 AliHLTITSVertexerZ::AliHLTITSVertexerZ():AliITSVertexerZ(){
35   // Constructor in case that there is no runloader
36
37   SetDiffPhiMax();
38   fX0 = 0;
39   fY0 = 0;
40   SetFirstLayerModules();
41   SetSecondLayerModules();
42   fZFound = 0;
43   fZsig = 0.;
44   //fITS = 0;
45   fZCombc = 0;
46   fZCombf = 0;
47   SetLowLimit();
48   SetHighLimit();
49   SetBinWidthCoarse();
50   SetBinWidthFine();
51   SetTolerance();
52   SetDebug();
53 }
54
55 AliHLTITSVertexerZ::AliHLTITSVertexerZ(TString filename,Float_t x0, Float_t y0):AliITSVertexerZ(filename,x0,y0)
56 {
57   // Standard Constructor
58 }
59
60 //______________________________________________________________________
61 AliESDVertex* AliHLTITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
62   // Defines the AliESDVertex for the current event
63
64   fCurrentVertex = 0;
65   AliRunLoader *rl =AliRunLoader::GetRunLoader();
66   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
67   itsLoader->LoadRecPoints();
68   rl->GetEvent(evnumber);
69
70   /*  if(!fITS)  {
71     fITS = (AliITS*)gAlice->GetModule("ITS");
72     if(!fITS) {
73       Error("FindVertexForCurrentEvent","AliITS object was not found");
74       return fCurrentVertex;
75     }
76   }
77   */
78
79   // fITS->SetTreeAddress();
80   rl->CdGAFile();
81
82   // AliITSgeom *geom = fITS->GetITSgeom();
83   AliITSgeom *geom = (AliITSgeom*)gDirectory->Get("AliITSgeom");
84   TTree *tR = itsLoader->TreeR();
85   
86   return FindVertexForCurrentEvent(geom,tR);
87 }
88
89 //______________________________________________________________________
90 AliESDVertex* AliHLTITSVertexerZ::FindVertexForCurrentEvent(AliITSgeom *geom,TTree *tR){
91   // Defines the AliESDVertex for the current event
92
93   fCurrentVertex = 0;
94
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.;
99
100   TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
101   TBranch *branch;
102   branch = tR->GetBranch("Clusters");
103   branch->SetAddress(&clusters);
104
105   Int_t nrpL1 = 0;
106   Int_t nrpL2 = 0;
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();
113   }
114   for(Int_t module= fFirstL2; module<=fLastL2;module++){
115     tR->GetEvent(module);
116     nrpL2+= clusters->GetEntriesFast();
117   }
118   if(nrpL1 == 0 || nrpL2 == 0){
119     ResetHistograms();
120     return fCurrentVertex;
121   }
122
123   Int_t nPhiBins = (Int_t)(TMath::Pi()/fDiffPhiMax);
124   Float_t phiBinSize = 2*TMath::Pi()/(Float_t)nPhiBins;
125
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];
143   }
144   
145   Float_t yshift = 0;
146   Float_t zshift[4] = {-10.708000, -3.536000, 3.536000, 10.708000};
147
148   yshift = 0.248499;
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);
159       gc[0]-=fX0;
160       gc[1]-=fY0;
161       Float_t xc1,yc1;
162       xc1=gc[0];
163       yc1=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);
167       if(bin>=nPhiBins || bin<0) bin = 0;
168       Int_t ind = ind1[bin];
169       if(ind<maxind1) {
170         phi1[bin][ind] = phi;
171         zc1[bin][ind]=gc[2]/fStepFine;
172         r1[bin][ind]=sqrt(xc1*xc1+yc1*yc1);
173         ind1[bin]++;
174       }
175     }
176     clusters->Delete();
177   }
178   yshift = 3.096207;
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);
188       gc[0]-=fX0;
189       gc[1]-=fY0;
190       Float_t xc2,yc2;
191       xc2=gc[0];
192       yc2=gc[1];
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];
198       if(ind<maxind2) {
199         phi2[bin][ind] = phi;
200         zc2[bin][ind]=gc[2]/fStepFine;
201         r2[bin][ind]=sqrt(xc2*xc2+yc2*yc2);
202         ind2[bin]++;
203       }
204     }
205     clusters->Delete();
206   }
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++) {
217         Int_t jbin2 = 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){
229               harray[bin]++;
230             }
231           }
232         }
233       }
234     }
235   }
236   for(Int_t i=0;i<nPhiBins;i++) {
237     delete [] zc1[i];
238     delete [] phi1[i];
239     delete [] r1[i];
240     delete [] zc2[i];
241     delete [] phi2[i];
242     delete [] r2[i];
243   }
244   delete [] zc1;
245   delete [] phi1;
246   delete [] r1;
247   delete [] ind1;
248   delete [] zc2;
249   delete [] phi2;
250   delete [] r2;
251   delete [] ind2;
252
253   Int_t nbinratio = (Int_t)(fStepCoarse/fStepFine+0.5);
254   Int_t nbincoarse = nbinfine/nbinratio;
255
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);
260   Stat_t contents=0;
261   Int_t counter=0;
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];
266     counter++;
267     if(counter==nbinratio) {
268       Int_t binc=bin/nbinratio; 
269       fZCombc->SetBinContent(binc+1,contents);
270       fZCombc->SetBinError(binc+1,TMath::Sqrt(contents));
271       contents=0;
272       counter=0;
273     }
274   }
275
276   delete [] harray;
277
278   if(fZCombc->GetEntries()==0){
279     Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
280     ResetHistograms();
281     return fCurrentVertex;
282   }
283   //  else {
284   //    cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
285   //  }
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));
290   Int_t niter = 0;
291   Bool_t goon = kTRUE;
292   Int_t num;
293   while(goon){
294     fZFound = 0.;
295     fZsig = 0.;
296     num=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);
301     }
302     if(num<2){
303       fZsig = 0.;
304     }
305     else {
306       Float_t radi =  fZsig/(num-1)-fZFound*fZFound/num/(num-1);
307       if(radi>0.)fZsig=TMath::Sqrt(radi);
308       else fZsig=0.;
309       fZFound/=num;
310     }
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));
314     niter++;
315     if(niter>=10){
316       goon = kFALSE;
317       Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
318     }
319   }
320   //  cout<<"Numer of Iterations "<<niter<<endl<<endl;
321   fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
322   fCurrentVertex->SetTitle("vertexer: HLT");
323   ResetHistograms();
324   PrintStatus();
325   return fCurrentVertex;
326 }