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