Moving the functions used to initialize TF1 and TF2 to the pivate part of the class
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerZ.cxx
CommitLineData
0c6af5c9 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
37ClassImp(AliITSVertexerZ)
38
39
40
41//______________________________________________________________________
42AliITSVertexerZ::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//______________________________________________________________________
62AliITSVertexerZ::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//______________________________________________________________________
83AliITSVertexerZ::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//______________________________________________________________________
90AliITSVertexerZ& 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//______________________________________________________________________
99AliITSVertexerZ::~AliITSVertexerZ() {
100 // Default Destructor
101 fITS = 0;
102 if(fZCombc)delete fZCombc;
103 if(fZCombf)delete fZCombf;
104}
105
106//______________________________________________________________________
d681bb2d 107AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
108 // Defines the AliESDVertex for the current event
0c6af5c9 109
110 fCurrentVertex = 0;
111 AliRunLoader *rl =AliRunLoader::GetRunLoader();
112 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
2257f27e 113 itsLoader->LoadRecPoints();
0c6af5c9 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
2257f27e 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 }
0c6af5c9 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;
2257f27e 159 // cout<<"Procesing module "<<module<<" ";
0c6af5c9 160 branch->GetEvent(module);
2257f27e 161 // cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
162 if(fUseV2Clusters){
163 Clusters2RecPoints(clusters,module,itsRec);
164 }
0c6af5c9 165 nrpL1+= itsRec->GetEntries();
166 fITS->ResetRecPoints();
167 }
168 for(Int_t module= fFirstL2; module<=fLastL2;module++){
169 branch->GetEvent(module);
2257f27e 170 if(fUseV2Clusters){
171 Clusters2RecPoints(clusters,module,itsRec);
172 }
0c6af5c9 173 nrpL2+= itsRec->GetEntries();
174 fITS->ResetRecPoints();
175 }
0c6af5c9 176 if(nrpL1 == 0 || nrpL2 == 0){
0c6af5c9 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);
2257f27e 192 if(fUseV2Clusters){
193 Clusters2RecPoints(clusters,module,itsRec);
194 }
0c6af5c9 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);
2257f27e 215 if(fUseV2Clusters){
216 Clusters2RecPoints(clusters,module,itsRec);
217 }
0c6af5c9 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 }
2257f27e 261 // else {
262 // cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
263 // }
0c6af5c9 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){
0c6af5c9 281 fZsig = 0.;
282 }
283 else {
979e3647 284 Float_t radi = fZsig/(num-1)-fZFound*fZFound/num/(num-1);
285 if(radi>0.)fZsig=TMath::Sqrt(radi);
286 else fZsig=0.;
0c6af5c9 287 fZFound/=num;
288 }
289 goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
0c6af5c9 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 }
2257f27e 298 // cout<<"Numer of Iterations "<<niter<<endl<<endl;
d681bb2d 299 fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
0c6af5c9 300 fCurrentVertex->SetTitle("vertexer: B");
301 ResetHistograms();
302 return fCurrentVertex;
303}
304
305//_____________________________________________________________________
306void 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//______________________________________________________________________
315void 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//________________________________________________________
337void 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