]>
Commit | Line | Data |
---|---|---|
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> | |
0c6af5c9 | 16 | #include <TString.h> |
0c6af5c9 | 17 | #include "AliITSLoader.h" |
0c6af5c9 | 18 | #include<TBranch.h> |
19 | #include<TClonesArray.h> | |
20 | #include<TH1.h> | |
0c6af5c9 | 21 | #include<TTree.h> |
22 | #include <AliITSgeom.h> | |
7d62fb64 | 23 | #include "AliITSDetTypeRec.h" |
0c6af5c9 | 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.; | |
0c6af5c9 | 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.; | |
0c6af5c9 | 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 | |
7d62fb64 | 96 | //fITS = 0; |
0c6af5c9 | 97 | if(fZCombc)delete fZCombc; |
98 | if(fZCombf)delete fZCombf; | |
99 | } | |
100 | ||
101 | //______________________________________________________________________ | |
d681bb2d | 102 | AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){ |
103 | // Defines the AliESDVertex for the current event | |
7d62fb64 | 104 | |
0c6af5c9 | 105 | fCurrentVertex = 0; |
106 | AliRunLoader *rl =AliRunLoader::GetRunLoader(); | |
7d62fb64 | 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 | ||
2257f27e | 113 | itsLoader->LoadRecPoints(); |
0c6af5c9 | 114 | rl->GetEvent(evnumber); |
115 | ||
7d62fb64 | 116 | AliITSDetTypeRec detTypeRec; |
0c6af5c9 | 117 | |
118 | TTree *tR = itsLoader->TreeR(); | |
7d62fb64 | 119 | detTypeRec.SetTreeAddressR(tR); |
0c6af5c9 | 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 | ||
7d62fb64 | 126 | itsRec = detTypeRec.RecPoints(); |
2257f27e | 127 | TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy; |
128 | TBranch *branch; | |
129 | if(fUseV2Clusters){ | |
130 | branch = tR->GetBranch("Clusters"); | |
131 | branch->SetAddress(&clusters); | |
132 | } | |
133 | else { | |
134 | branch = tR->GetBranch("ITSRecPoints"); | |
135 | } | |
0c6af5c9 | 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; | |
2257f27e | 148 | // cout<<"Procesing module "<<module<<" "; |
0c6af5c9 | 149 | branch->GetEvent(module); |
2257f27e | 150 | // cout<<"Number of clusters "<<clusters->GetEntries()<<endl; |
151 | if(fUseV2Clusters){ | |
152 | Clusters2RecPoints(clusters,module,itsRec); | |
153 | } | |
0c6af5c9 | 154 | nrpL1+= itsRec->GetEntries(); |
7d62fb64 | 155 | detTypeRec.ResetRecPoints(); |
0c6af5c9 | 156 | } |
157 | for(Int_t module= fFirstL2; module<=fLastL2;module++){ | |
158 | branch->GetEvent(module); | |
2257f27e | 159 | if(fUseV2Clusters){ |
160 | Clusters2RecPoints(clusters,module,itsRec); | |
161 | } | |
0c6af5c9 | 162 | nrpL2+= itsRec->GetEntries(); |
7d62fb64 | 163 | detTypeRec.ResetRecPoints(); |
0c6af5c9 | 164 | } |
0c6af5c9 | 165 | if(nrpL1 == 0 || nrpL2 == 0){ |
0c6af5c9 | 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); | |
2257f27e | 181 | if(fUseV2Clusters){ |
182 | Clusters2RecPoints(clusters,module,itsRec); | |
183 | } | |
0c6af5c9 | 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->GetX(); | |
188 | lc[2]=recp->GetZ(); | |
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 | } | |
7d62fb64 | 199 | detTypeRec.ResetRecPoints(); |
0c6af5c9 | 200 | } |
201 | ind = 0; | |
202 | for(Int_t module= fFirstL2; module<=fLastL2;module++){ | |
203 | branch->GetEvent(module); | |
2257f27e | 204 | if(fUseV2Clusters){ |
205 | Clusters2RecPoints(clusters,module,itsRec); | |
206 | } | |
0c6af5c9 | 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->GetX(); | |
211 | lc[2]=recp->GetZ(); | |
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 | } | |
7d62fb64 | 222 | detTypeRec.ResetRecPoints(); |
0c6af5c9 | 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 | } | |
2257f27e | 250 | // else { |
251 | // cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl; | |
252 | // } | |
0c6af5c9 | 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){ | |
0c6af5c9 | 270 | fZsig = 0.; |
271 | } | |
272 | else { | |
979e3647 | 273 | Float_t radi = fZsig/(num-1)-fZFound*fZFound/num/(num-1); |
274 | if(radi>0.)fZsig=TMath::Sqrt(radi); | |
275 | else fZsig=0.; | |
0c6af5c9 | 276 | fZFound/=num; |
277 | } | |
278 | goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance; | |
0c6af5c9 | 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 | } | |
2257f27e | 287 | // cout<<"Numer of Iterations "<<niter<<endl<<endl; |
d681bb2d | 288 | fCurrentVertex = new AliESDVertex(fZFound,fZsig,num); |
0c6af5c9 | 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 |