Moving ALiITSVertex to AliESDVertex
[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->ReloadRecPoints();
114
115   rl->GetEvent(evnumber);
116
117   if(!fITS)  {
118     fITS = (AliITS*)gAlice->GetModule("ITS");
119     if(!fITS) {
120       Error("FindVertexForCurrentEvent","AliITS object was not found");
121       return fCurrentVertex;
122     }
123   }
124
125   fITS->SetTreeAddress();
126
127   AliITSgeom *geom = fITS->GetITSgeom();
128
129   TTree *tR = itsLoader->TreeR();
130   TClonesArray *itsRec  = 0;
131   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
132   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
133   Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
134   Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
135
136   itsRec  = fITS->RecPoints();
137   TBranch *branch = tR->GetBranch("ITSRecPoints");
138
139   Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
140   Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
141   if(fZCombc)delete fZCombc;
142   fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
143   if(fZCombf)delete fZCombf;
144   fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
145
146   Int_t nrpL1 = 0;
147   Int_t nrpL2 = 0;
148   for(Int_t module= fFirstL1; module<=fLastL1;module++){
149     if(module%4==0 || module%4==3)continue;
150     branch->GetEvent(module);
151     nrpL1+= itsRec->GetEntries();
152     fITS->ResetRecPoints();
153   }
154   for(Int_t module= fFirstL2; module<=fLastL2;module++){
155     branch->GetEvent(module);
156     nrpL2+= itsRec->GetEntries();
157     fITS->ResetRecPoints();
158   }
159   cout<<"nrpL1 = "<<nrpL1<<"; nrpL2= "<<nrpL2<<endl;
160   if(nrpL1 == 0 || nrpL2 == 0){
161     cout<<"achi, achi\n";
162     ResetHistograms();
163     return fCurrentVertex;
164   }
165   Float_t *xc1 = new Float_t [nrpL1];
166   Float_t *yc1 = new Float_t [nrpL1];
167   Float_t *zc1 = new Float_t [nrpL1];
168   Float_t *phi1 = new Float_t [nrpL1];
169   Float_t *xc2 = new Float_t [nrpL2];
170   Float_t *yc2 = new Float_t [nrpL2];
171   Float_t *zc2 = new Float_t [nrpL2];
172   Float_t *phi2 = new Float_t [nrpL2];
173   Int_t ind = 0;
174   for(Int_t module= fFirstL1; module<=fLastL1;module++){
175     if(module%4==0 || module%4==3)continue;
176     branch->GetEvent(module);
177     Int_t nrecp1 = itsRec->GetEntries();
178     for(Int_t j=0;j<nrecp1;j++){
179       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
180       lc[0]=recp->GetX();
181       lc[2]=recp->GetZ();
182       geom->LtoG(module,lc,gc);
183       gc[0]-=fX0;
184       gc[1]-=fY0;
185       xc1[ind]=gc[0];
186       yc1[ind]=gc[1];
187       zc1[ind]=gc[2];
188       phi1[ind] = TMath::ATan2(gc[1],gc[0]);
189       if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
190       ind++;
191     }
192     fITS->ResetRecPoints();
193   }
194   ind = 0;
195   for(Int_t module= fFirstL2; module<=fLastL2;module++){
196     branch->GetEvent(module);
197     Int_t nrecp2 = itsRec->GetEntries();
198     for(Int_t j=0;j<nrecp2;j++){
199       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
200       lc[0]=recp->GetX();
201       lc[2]=recp->GetZ();
202       geom->LtoG(module,lc,gc);
203       gc[0]-=fX0;
204       gc[1]-=fY0;
205       xc2[ind]=gc[0];
206       yc2[ind]=gc[1];
207       zc2[ind]=gc[2];
208       phi2[ind] = TMath::ATan2(gc[1],gc[0]);
209       if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
210       ind++;
211     }
212     fITS->ResetRecPoints();
213   }
214   for(Int_t i=0;i<nrpL1;i++){
215     Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]);
216     for(Int_t j=0;j<nrpL2;j++){
217       Float_t diff = TMath::Abs(phi2[j]-phi1[i]);
218       if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
219       if(diff<fDiffPhiMax){
220         Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]);
221         Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1);
222         fZCombf->Fill(zr0);
223         fZCombc->Fill(zr0);
224       }
225     }
226   }
227   delete [] xc1;
228   delete [] yc1;
229   delete [] zc1;
230   delete [] phi1;
231   delete [] xc2;
232   delete [] yc2;
233   delete [] zc2;
234   delete [] phi2;
235   if(fZCombc->GetEntries()==0){
236     Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
237     ResetHistograms();
238     return fCurrentVertex;
239   }
240   else {
241     cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
242   }
243   Int_t bi = fZCombc->GetMaximumBin();
244   Float_t centre = fZCombc->GetBinCenter(bi);
245   Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
246   Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
247   Int_t niter = 0;
248   Bool_t goon = kTRUE;
249   Int_t num;
250   while(goon){
251     fZFound = 0.;
252     fZsig = 0.;
253     num=0;
254     for(Int_t n=n1;n<=n2;n++){
255       fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
256       num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
257       fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
258     }
259     if(num<2){
260       fZsig = 0.;
261     }
262     else {
263       Float_t radi =  fZsig/(num-1)-fZFound*fZFound/num/(num-1);
264       if(radi>0.)fZsig=TMath::Sqrt(radi);
265       else fZsig=0.;
266       fZFound/=num;
267     }
268     goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
269     cout<<"Simmetria "<<TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))<<endl;
270     n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
271     n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
272     niter++;
273     if(niter>=10){
274       goon = kFALSE;
275       Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
276     }
277   }
278   cout<<"Numer of Iterations "<<niter<<endl<<endl;
279   fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
280   fCurrentVertex->SetTitle("vertexer: B");
281   ResetHistograms();
282   return fCurrentVertex;
283 }
284
285 //_____________________________________________________________________
286 void AliITSVertexerZ::ResetHistograms(){
287   // delete TH1 data members
288   if(fZCombc)delete fZCombc;
289   if(fZCombf)delete fZCombf;
290   fZCombc = 0;
291   fZCombf = 0;
292 }
293
294 //______________________________________________________________________
295 void AliITSVertexerZ::FindVertices(){
296   // computes the vertices of the events in the range FirstEvent - LastEvent
297   AliRunLoader *rl = AliRunLoader::GetRunLoader();
298   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
299   itsLoader->ReloadRecPoints();
300   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
301     cout<<"Processing event "<<i<<endl;
302     rl->GetEvent(i);
303     FindVertexForCurrentEvent(i);
304     if(fCurrentVertex){
305       WriteCurrentVertex();
306     }
307     else {
308       if(fDebug>0){
309         cout<<"Vertex not found for event "<<i<<endl;
310         cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
311       }
312     }
313   }
314 }
315
316 //________________________________________________________
317 void AliITSVertexerZ::PrintStatus() const {
318   // Print current status
319   cout <<"=======================================================\n";
320   cout <<" First layer first and last modules: "<<fFirstL1<<", ";
321   cout <<fLastL1<<endl;
322   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
323   cout <<fLastL2<<endl;
324   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
325   cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
326   cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
327   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
328   cout <<" Debug flag: "<<fDebug<<endl;
329   cout <<"First event to be processed "<<fFirstEvent;
330   cout <<"\n Last event to be processed "<<fLastEvent<<endl;
331  
332  
333   cout <<"=======================================================\n";
334 }
335