Primary vertex reconstruction and standalone ITS tracking in the reconstruction chain
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerZ.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2003, 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 <AliITSVertexerZ.h>
16 #include <Riostream.h>
17 #include <TString.h>
18 #include <AliITS.h>
19 #include "AliITSLoader.h"
20 #include <AliRun.h>
21 #include<TBranch.h>
22 #include<TClonesArray.h>
23 #include<TH1.h>
24 #include<TMath.h>
25 #include<TTree.h>
26 #include <AliITSgeom.h>
27 #include <AliITSRecPoint.h>
28
29 /////////////////////////////////////////////////////////////////
30 // this class implements a fast method to determine
31 // the Z coordinate of the primary vertex
32 // for p-p collisions it seems to give comparable or better results
33 // with respect to what obtained with AliITSVertexerPPZ
34 // It can be used successfully with Pb-Pb collisions
35 ////////////////////////////////////////////////////////////////
36
37 ClassImp(AliITSVertexerZ)
38
39
40
41 //______________________________________________________________________
42 AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer() {
43   // Default constructor
44   SetDiffPhiMax(0);
45   fX0 = 0.;
46   fY0 = 0.;
47   SetFirstLayerModules(0);
48   SetSecondLayerModules(0);
49   fZFound = 0;
50   fZsig = 0.;
51   fITS = 0;
52   fZCombc = 0;
53   fZCombf = 0;
54   SetLowLimit(0.);
55   SetHighLimit(0.);
56   SetBinWidthCoarse(0.);
57   SetBinWidthFine(0.);
58   SetTolerance(0.);
59 }
60
61 //______________________________________________________________________
62 AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
63   // Standard Constructor
64   SetDiffPhiMax();
65   fX0 = x0;
66   fY0 = y0;
67   SetFirstLayerModules();
68   SetSecondLayerModules();
69   fZFound = 0;
70   fZsig = 0.;
71   fITS = 0;
72   fZCombc = 0;
73   fZCombf = 0;
74   SetLowLimit();
75   SetHighLimit();
76   SetBinWidthCoarse();
77   SetBinWidthFine();
78   SetTolerance();
79
80 }
81
82 //______________________________________________________________________
83 AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr) {
84   // Copy constructor
85   // Copies are not allowed. The method is protected to avoid misuse.
86   Error("AliITSVertexerZ","Copy constructor not allowed\n");
87 }
88
89 //______________________________________________________________________
90 AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ& /* vtxr */){
91   // Assignment operator
92   // Assignment is not allowed. The method is protected to avoid misuse.
93   Error("= operator","Assignment operator not allowed\n");
94   return *this;
95 }
96
97
98 //______________________________________________________________________
99 AliITSVertexerZ::~AliITSVertexerZ() {
100   // Default Destructor
101   fITS = 0;
102   if(fZCombc)delete fZCombc;
103   if(fZCombf)delete fZCombf;
104 }
105
106 //______________________________________________________________________
107 AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
108   // Defines the AliESDVertex for the current event
109
110   fCurrentVertex = 0;
111   AliRunLoader *rl =AliRunLoader::GetRunLoader();
112   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
113   itsLoader->LoadRecPoints();
114   rl->GetEvent(evnumber);
115
116   if(!fITS)  {
117     fITS = (AliITS*)gAlice->GetModule("ITS");
118     if(!fITS) {
119       Error("FindVertexForCurrentEvent","AliITS object was not found");
120       return fCurrentVertex;
121     }
122   }
123
124   fITS->SetTreeAddress();
125
126   AliITSgeom *geom = fITS->GetITSgeom();
127
128   TTree *tR = itsLoader->TreeR();
129   TClonesArray *itsRec  = 0;
130   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
131   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
132   Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
133   Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
134
135   itsRec = fITS->RecPoints();
136
137   //cout<<"Address of itsRec = "<<itsRec<<endl;
138   TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
139   TBranch *branch;
140   if(fUseV2Clusters){
141     branch = tR->GetBranch("Clusters");
142     branch->SetAddress(&clusters);
143   }
144   else {
145     branch = tR->GetBranch("ITSRecPoints");
146   }
147
148   Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
149   Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
150   if(fZCombc)delete fZCombc;
151   fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
152   if(fZCombf)delete fZCombf;
153   fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
154
155   Int_t nrpL1 = 0;
156   Int_t nrpL2 = 0;
157   for(Int_t module= fFirstL1; module<=fLastL1;module++){
158     if(module%4==0 || module%4==3)continue;
159     //   cout<<"Procesing module "<<module<<" ";
160     branch->GetEvent(module);
161     //    cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
162     if(fUseV2Clusters){
163       Clusters2RecPoints(clusters,module,itsRec);
164     }
165     nrpL1+= itsRec->GetEntries();
166     fITS->ResetRecPoints();
167   }
168   for(Int_t module= fFirstL2; module<=fLastL2;module++){
169     branch->GetEvent(module);
170     if(fUseV2Clusters){
171       Clusters2RecPoints(clusters,module,itsRec);
172     }
173     nrpL2+= itsRec->GetEntries();
174     fITS->ResetRecPoints();
175   }
176   if(nrpL1 == 0 || nrpL2 == 0){
177     ResetHistograms();
178     return fCurrentVertex;
179   }
180   Float_t *xc1 = new Float_t [nrpL1];
181   Float_t *yc1 = new Float_t [nrpL1];
182   Float_t *zc1 = new Float_t [nrpL1];
183   Float_t *phi1 = new Float_t [nrpL1];
184   Float_t *xc2 = new Float_t [nrpL2];
185   Float_t *yc2 = new Float_t [nrpL2];
186   Float_t *zc2 = new Float_t [nrpL2];
187   Float_t *phi2 = new Float_t [nrpL2];
188   Int_t ind = 0;
189   for(Int_t module= fFirstL1; module<=fLastL1;module++){
190     if(module%4==0 || module%4==3)continue;
191     branch->GetEvent(module);
192     if(fUseV2Clusters){
193       Clusters2RecPoints(clusters,module,itsRec);
194     }
195     Int_t nrecp1 = itsRec->GetEntries();
196     for(Int_t j=0;j<nrecp1;j++){
197       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
198       lc[0]=recp->GetX();
199       lc[2]=recp->GetZ();
200       geom->LtoG(module,lc,gc);
201       gc[0]-=fX0;
202       gc[1]-=fY0;
203       xc1[ind]=gc[0];
204       yc1[ind]=gc[1];
205       zc1[ind]=gc[2];
206       phi1[ind] = TMath::ATan2(gc[1],gc[0]);
207       if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
208       ind++;
209     }
210     fITS->ResetRecPoints();
211   }
212   ind = 0;
213   for(Int_t module= fFirstL2; module<=fLastL2;module++){
214     branch->GetEvent(module);
215     if(fUseV2Clusters){
216       Clusters2RecPoints(clusters,module,itsRec);
217     }
218     Int_t nrecp2 = itsRec->GetEntries();
219     for(Int_t j=0;j<nrecp2;j++){
220       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
221       lc[0]=recp->GetX();
222       lc[2]=recp->GetZ();
223       geom->LtoG(module,lc,gc);
224       gc[0]-=fX0;
225       gc[1]-=fY0;
226       xc2[ind]=gc[0];
227       yc2[ind]=gc[1];
228       zc2[ind]=gc[2];
229       phi2[ind] = TMath::ATan2(gc[1],gc[0]);
230       if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
231       ind++;
232     }
233     fITS->ResetRecPoints();
234   }
235   for(Int_t i=0;i<nrpL1;i++){
236     Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]);
237     for(Int_t j=0;j<nrpL2;j++){
238       Float_t diff = TMath::Abs(phi2[j]-phi1[i]);
239       if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
240       if(diff<fDiffPhiMax){
241         Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]);
242         Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1);
243         fZCombf->Fill(zr0);
244         fZCombc->Fill(zr0);
245       }
246     }
247   }
248   delete [] xc1;
249   delete [] yc1;
250   delete [] zc1;
251   delete [] phi1;
252   delete [] xc2;
253   delete [] yc2;
254   delete [] zc2;
255   delete [] phi2;
256   if(fZCombc->GetEntries()==0){
257     Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
258     ResetHistograms();
259     return fCurrentVertex;
260   }
261   //  else {
262   //    cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
263   //  }
264   Int_t bi = fZCombc->GetMaximumBin();
265   Float_t centre = fZCombc->GetBinCenter(bi);
266   Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
267   Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
268   Int_t niter = 0;
269   Bool_t goon = kTRUE;
270   Int_t num;
271   while(goon){
272     fZFound = 0.;
273     fZsig = 0.;
274     num=0;
275     for(Int_t n=n1;n<=n2;n++){
276       fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
277       num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
278       fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
279     }
280     if(num<2){
281       fZsig = 0.;
282     }
283     else {
284       Float_t radi =  fZsig/(num-1)-fZFound*fZFound/num/(num-1);
285       if(radi>0.)fZsig=TMath::Sqrt(radi);
286       else fZsig=0.;
287       fZFound/=num;
288     }
289     goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
290     n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
291     n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
292     niter++;
293     if(niter>=10){
294       goon = kFALSE;
295       Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
296     }
297   }
298   //  cout<<"Numer of Iterations "<<niter<<endl<<endl;
299   fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
300   fCurrentVertex->SetTitle("vertexer: B");
301   ResetHistograms();
302   return fCurrentVertex;
303 }
304
305 //_____________________________________________________________________
306 void AliITSVertexerZ::ResetHistograms(){
307   // delete TH1 data members
308   if(fZCombc)delete fZCombc;
309   if(fZCombf)delete fZCombf;
310   fZCombc = 0;
311   fZCombf = 0;
312 }
313
314 //______________________________________________________________________
315 void AliITSVertexerZ::FindVertices(){
316   // computes the vertices of the events in the range FirstEvent - LastEvent
317   AliRunLoader *rl = AliRunLoader::GetRunLoader();
318   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
319   itsLoader->ReloadRecPoints();
320   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
321     cout<<"Processing event "<<i<<endl;
322     rl->GetEvent(i);
323     FindVertexForCurrentEvent(i);
324     if(fCurrentVertex){
325       WriteCurrentVertex();
326     }
327     else {
328       if(fDebug>0){
329         cout<<"Vertex not found for event "<<i<<endl;
330         cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
331       }
332     }
333   }
334 }
335
336 //________________________________________________________
337 void AliITSVertexerZ::PrintStatus() const {
338   // Print current status
339   cout <<"=======================================================\n";
340   cout <<" First layer first and last modules: "<<fFirstL1<<", ";
341   cout <<fLastL1<<endl;
342   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
343   cout <<fLastL2<<endl;
344   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
345   cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
346   cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
347   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
348   cout <<" Debug flag: "<<fDebug<<endl;
349   cout <<"First event to be processed "<<fFirstEvent;
350   cout <<"\n Last event to be processed "<<fLastEvent<<endl;
351  
352  
353   cout <<"=======================================================\n";
354 }
355