]>
Commit | Line | Data |
---|---|---|
f4f84117 | 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 "AliITSVertexerZD.h" | |
16 | #include<TBranch.h> | |
17 | #include<TClonesArray.h> | |
18 | #include<TH1.h> | |
19 | #include <TString.h> | |
20 | #include<TTree.h> | |
21 | #include "AliESDVertex.h" | |
22 | #include "AliITSgeomTGeo.h" | |
23 | #include "AliITSRecPoint.h" | |
24 | #include "AliITSRecPointContainer.h" | |
25 | #include "AliITSZPoint.h" | |
26 | ||
27 | ///////////////////////////////////////////////////////////////// | |
28 | // this class implements a fast method to determine | |
29 | // the Z coordinate of the primary vertex | |
30 | // it is based on the 2 layers of SDD detectors | |
31 | // wider Z acceptance w.r.t. AliITSVertexerZ | |
32 | // lower resolution | |
33 | // origin masera@to.infn.it | |
34 | //////////////////////////////////////////////////////////////// | |
35 | ||
36 | using std::endl; | |
37 | using std::cout; | |
38 | ClassImp(AliITSVertexerZD) | |
39 | ||
40 | ||
41 | ||
42 | //______________________________________________________________________ | |
43 | AliITSVertexerZD::AliITSVertexerZD():AliITSVertexer(), | |
44 | fFirstL1(0), | |
45 | fLastL1(0), | |
46 | fFirstL2(0), | |
47 | fLastL2(0), | |
48 | fDiffPhiMax(0), | |
49 | fZFound(0), | |
50 | fZsig(0.), | |
51 | fZCombc(0), | |
52 | fLowLim(0.), | |
53 | fHighLim(0.), | |
54 | fStepCoarse(0), | |
55 | fTolerance(0.), | |
56 | fMaxIter(0), | |
57 | fWindowWidth(0), | |
58 | fSearchForPileup(kTRUE) | |
59 | { | |
60 | // Default constructor | |
61 | SetDiffPhiMax(); | |
62 | SetFirstLayerModules(); | |
63 | SetSecondLayerModules(); | |
64 | SetLowLimit(); | |
65 | SetHighLimit(); | |
66 | SetBinWidthCoarse(); | |
67 | SetTolerance(); | |
68 | SetPPsetting(); | |
69 | ConfigIterations(); | |
70 | SetWindowWidth(); | |
71 | } | |
72 | ||
73 | //______________________________________________________________________ | |
74 | AliITSVertexerZD::AliITSVertexerZD(Float_t x0, Float_t y0):AliITSVertexer(), | |
75 | fFirstL1(0), | |
76 | fLastL1(0), | |
77 | fFirstL2(0), | |
78 | fLastL2(0), | |
79 | fDiffPhiMax(0), | |
80 | fZFound(0), | |
81 | fZsig(0.), | |
82 | fZCombc(0), | |
83 | fLowLim(0.), | |
84 | fHighLim(0.), | |
85 | fStepCoarse(0), | |
86 | fTolerance(0.), | |
87 | fMaxIter(0), | |
88 | fWindowWidth(0), | |
89 | fSearchForPileup(kTRUE) | |
90 | { | |
91 | // Standard Constructor | |
92 | SetDiffPhiMax(); | |
93 | SetFirstLayerModules(); | |
94 | SetSecondLayerModules(); | |
95 | SetLowLimit(); | |
96 | SetHighLimit(); | |
97 | SetBinWidthCoarse(); | |
98 | SetTolerance(); | |
99 | SetPPsetting(); | |
100 | ConfigIterations(); | |
101 | SetWindowWidth(); | |
102 | SetVtxStart((Double_t)x0,(Double_t)y0,0.); | |
103 | ||
104 | } | |
105 | ||
106 | //______________________________________________________________________ | |
107 | AliITSVertexerZD::~AliITSVertexerZD() { | |
108 | // Destructor | |
109 | delete fZCombc; | |
110 | } | |
111 | ||
112 | //______________________________________________________________________ | |
113 | void AliITSVertexerZD::ConfigIterations(Int_t noiter,Float_t *ptr){ | |
114 | // configure the iterative procedure to gain efficiency for | |
115 | // pp events with very low multiplicity | |
116 | Float_t defaults[5]={0.02,0.05,0.1,0.2,0.3}; | |
117 | fMaxIter=noiter; | |
118 | if(noiter>5){ | |
119 | Error("ConfigIterations","Maximum number of iterations is 5\n"); | |
120 | fMaxIter=5; | |
121 | } | |
122 | for(Int_t j=0;j<5;j++)fPhiDiffIter[j]=defaults[j]; | |
123 | if(ptr)for(Int_t j=0;j<fMaxIter;j++)fPhiDiffIter[j]=ptr[j]; | |
124 | } | |
125 | ||
126 | //______________________________________________________________________ | |
127 | Int_t AliITSVertexerZD::GetPeakRegion(TH1F*h, Int_t &binmin, Int_t &binmax){ | |
128 | // Finds a region around a peak in the Z histogram | |
129 | // Case of 2 peaks is treated | |
130 | Int_t imax=h->GetNbinsX(); | |
131 | Float_t maxval=0; | |
132 | Int_t bi1=h->GetMaximumBin(); | |
133 | Int_t bi2=0; | |
134 | for(Int_t i=imax;i>=1;i--){ | |
135 | if(h->GetBinContent(i)>maxval){ | |
136 | maxval=h->GetBinContent(i); | |
137 | bi2=i; | |
138 | } | |
139 | } | |
140 | Int_t npeaks=0; | |
141 | ||
142 | if(bi1==bi2){ | |
143 | binmin=bi1-3; | |
144 | binmax=bi1+3; | |
145 | npeaks=1; | |
146 | }else{ | |
147 | TH1F *copy = new TH1F(*h); | |
148 | copy->SetBinContent(bi1,0.); | |
149 | copy->SetBinContent(bi2,0.); | |
150 | Int_t l1=TMath::Max(bi1-3,1); | |
151 | Int_t l2=TMath::Min(bi1+3,h->GetNbinsX()); | |
152 | Float_t cont1=copy->Integral(l1,l2); | |
153 | Int_t ll1=TMath::Max(bi2-3,1); | |
154 | Int_t ll2=TMath::Min(bi2+3,h->GetNbinsX()); | |
155 | Float_t cont2=copy->Integral(ll1,ll2); | |
156 | if(cont1>cont2){ | |
157 | binmin=l1; | |
158 | binmax=l2; | |
159 | npeaks=1; | |
160 | } | |
161 | if(cont2>cont1){ | |
162 | binmin=ll1; | |
163 | binmax=ll2; | |
164 | npeaks=1; | |
165 | } | |
166 | if(cont1==cont2){ | |
167 | binmin=l1; | |
168 | binmax=ll2; | |
169 | if(bi2-bi1==1) npeaks=1; | |
170 | else npeaks=2; | |
171 | } | |
172 | delete copy; | |
173 | } | |
174 | return npeaks; | |
175 | } | |
176 | //______________________________________________________________________ | |
177 | AliESDVertex* AliITSVertexerZD::FindVertexForCurrentEvent(TTree *itsClusterTree){ | |
178 | // Defines the AliESDVertex for the current event | |
179 | VertexZFinder(itsClusterTree); | |
180 | Int_t ntrackl=0; | |
181 | for(Int_t iteraz=0;iteraz<fMaxIter;iteraz++){ | |
182 | if(fCurrentVertex) ntrackl=fCurrentVertex->GetNContributors(); | |
183 | if(!fCurrentVertex || ntrackl==0 || ntrackl==-1){ | |
184 | Float_t diffPhiMaxOrig=fDiffPhiMax; | |
185 | fDiffPhiMax=GetPhiMaxIter(iteraz); | |
186 | VertexZFinder(itsClusterTree); | |
187 | fDiffPhiMax=diffPhiMaxOrig; | |
188 | } | |
189 | } | |
190 | if(fComputeMultiplicity) FindMultiplicity(itsClusterTree); | |
191 | return fCurrentVertex; | |
192 | } | |
193 | ||
194 | //______________________________________________________________________ | |
195 | void AliITSVertexerZD::VertexZFinder(TTree *itsClusterTree){ | |
196 | // Defines the AliESDVertex for the current event | |
197 | //debug printf("=========== VertexZFinder has been called \n"); | |
198 | fCurrentVertex = 0; | |
199 | Double_t startPos[3]={GetNominalPos()[0],GetNominalPos()[1],GetNominalPos()[2]}; | |
200 | Double_t startCov[6]={GetNominalCov()[0],GetNominalCov()[1],GetNominalCov()[2], | |
201 | GetNominalCov()[3],GetNominalCov()[4],GetNominalCov()[5]}; | |
202 | ResetVertex(); | |
203 | TClonesArray *itsRec = 0; | |
204 | // lc1 and gc1 are local and global coordinates for layer 1 | |
205 | Float_t gc1[3]={0.,0.,0.}; // ; for(Int_t ii=0; ii<3; ii++) gc1[ii]=0.; | |
206 | // lc2 and gc2 are local and global coordinates for layer 2 | |
207 | Float_t gc2[3]={0.,0.,0.}; //; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.; | |
208 | AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance(); | |
209 | rpcont->FetchClusters(0,itsClusterTree); | |
210 | if(!rpcont->IsSDDActive()){ | |
211 | AliWarning("Null pointer for RecPoints branch, vertex not calculated"); | |
212 | ResetHistograms(); | |
213 | fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-2); | |
214 | return; | |
215 | } | |
216 | ||
217 | Int_t nrpL1 = 0; | |
218 | Int_t nrpL2 = 0; | |
219 | nrpL1=rpcont->GetNClustersInLayerFast(3); | |
220 | nrpL2=rpcont->GetNClustersInLayerFast(4); | |
221 | ||
222 | if(nrpL1 == 0 || nrpL2 == 0){ | |
223 | AliDebug(1,Form("No RecPoints in at least one SDD layer (%d %d)",nrpL1,nrpL2)); | |
224 | ResetHistograms(); | |
225 | fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-2); | |
226 | return; | |
227 | } | |
228 | // Force a coarse bin size of 200 microns if the number of clusters on layer 2 | |
229 | // is low | |
230 | if(nrpL2<fPPsetting[0])SetBinWidthCoarse(fPPsetting[1]); | |
231 | // By default nbincoarse=(10+10)/0.01=2000 | |
232 | Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse); | |
233 | if(fZCombc)delete fZCombc; | |
234 | fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse); | |
235 | ||
236 | ||
237 | Int_t nEntriesMod[260]; | |
238 | TClonesArray* recpArr[260]; | |
239 | for(Int_t modul=0; modul<260; ++modul) { | |
240 | recpArr[modul]=rpcont->UncheckedGetClusters(modul+240); | |
241 | nEntriesMod[modul]=recpArr[modul]->GetEntriesFast(); | |
242 | } | |
243 | ||
244 | Int_t maxdim=TMath::Min(nrpL1*nrpL2,50000); // temporary; to limit the size in PbPb | |
245 | static TClonesArray points("AliITSZPoint",maxdim); | |
246 | Int_t nopoints =0; | |
247 | for(Int_t modul1= fFirstL1; modul1<=fLastL1;modul1++){ // Loop on modules of layer 1 | |
248 | // UShort_t ladder=int(modul1/6)+1; // ladders are numbered starting from 1 | |
249 | TClonesArray *prpl1=recpArr[modul1-240]; //rpcont->UncheckedGetClusters(modul1); | |
250 | Int_t nrecp1 = nEntriesMod[modul1-240]; //prpl1->GetEntries(); | |
251 | for(Int_t j1=0;j1<nrecp1;j1++){ | |
252 | AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1->At(j1); | |
253 | recp1->GetGlobalXYZ(gc1); | |
254 | gc1[0]-=GetNominalPos()[0]; // Possible beam offset in the bending plane | |
255 | gc1[1]-=GetNominalPos()[1]; // " " | |
256 | Float_t phi1 = TMath::ATan2(gc1[1],gc1[0]); | |
257 | if(phi1<0)phi1+=TMath::TwoPi(); | |
258 | for(Int_t modul2=fFirstL2;modul2<=fLastL2;modul2++){ | |
259 | itsRec=recpArr[modul2-240]; // rpcont->UncheckedGetClusters(modul2); | |
260 | Int_t nrecp2 = nEntriesMod[modul2-240]; // itsRec->GetEntries(); | |
261 | for(Int_t j2=0;j2<nrecp2;j2++){ | |
262 | AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2); | |
263 | recp2->GetGlobalXYZ(gc2); | |
264 | gc2[0]-=GetNominalPos()[0]; | |
265 | gc2[1]-=GetNominalPos()[1]; | |
266 | Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]); | |
267 | if(phi2<0)phi2+=TMath::TwoPi(); | |
268 | ||
269 | Float_t diff = TMath::Abs(phi2-phi1); | |
270 | if(diff>TMath::Pi())diff=TMath::TwoPi()-diff; | |
271 | if(diff<fDiffPhiMax){ | |
272 | //================= DEBUG | |
273 | /* Bool_t goodMatch = kFALSE; | |
274 | for(Int_t kl1=0;kl1<3;kl1++){ | |
275 | Int_t label1=recp1->GetLabel(kl1); | |
276 | if(label1>0){ | |
277 | for(Int_t kl2=0;kl2<3;kl2++){ | |
278 | if(label1 == recp2->GetLabel(kl2))goodMatch = kTRUE; | |
279 | } | |
280 | } | |
281 | } | |
282 | */ | |
283 | //================= END DEBUG | |
284 | Float_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); | |
285 | Float_t zc1=gc1[2]; | |
286 | Float_t erz1=recp1->GetSigmaZ2(); | |
287 | Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); | |
288 | //debug printf("r1 =%f ; r2=%f \n",r1,r2); | |
289 | Float_t zc2=gc2[2]; | |
290 | Float_t erz2=recp2->GetSigmaZ2(); | |
291 | // Float_t tgth=(zc2[j]-zc1[i])/(r2-r1); // slope (used for multiple scattering) | |
292 | Float_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius | |
293 | Float_t ezr0q=(r2*r2*erz1+r1*r1*erz2)/((r2-r1)*(r2-r1)); //error on Z @ null radius | |
294 | /* | |
295 | if(goodMatch)printf("Good tracklet: z0=%f, ezr0q=%f\n",zr0,ezr0q); | |
296 | else printf("Fake tracklet: z0=%f, ezr0q=%f\n",zr0,ezr0q); | |
297 | if(goodMatch){ | |
298 | printf("Labels: "); | |
299 | for(Int_t kkk=0;kkk<3;kkk++)printf("%d ",recp1->GetLabel(kkk)); | |
300 | for(Int_t kkk=0;kkk<3;kkk++)printf("%d ",recp1->GetLabel(kkk)); | |
301 | } | |
302 | printf("\n"); | |
303 | if(nopoints<maxdim) new(points[nopoints++])AliITSZPoint(zr0,ezr0q,goodMatch); | |
304 | */ | |
305 | if(nopoints<maxdim) new(points[nopoints++])AliITSZPoint(zr0,ezr0q); | |
306 | fZCombc->Fill(zr0); | |
307 | } | |
308 | } | |
309 | } | |
310 | } | |
311 | } | |
312 | points.Sort(); | |
313 | ||
314 | Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1); | |
315 | if(contents<1.){ | |
316 | // Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n"); | |
317 | ResetHistograms(); | |
318 | fCurrentVertex = new AliESDVertex(startPos,startCov,99999.,-1); | |
319 | points.Clear(); | |
320 | return; | |
321 | } | |
322 | ||
323 | TH1F *hc = fZCombc; | |
324 | ||
325 | ||
326 | if(hc->GetBinContent(hc->GetMaximumBin())<3)hc->Rebin(4); | |
327 | Int_t binmin,binmax; | |
328 | Int_t nPeaks=GetPeakRegion(hc,binmin,binmax); | |
329 | if(nPeaks==2)AliDebug(2,"2 peaks found"); | |
330 | Float_t zm =0.; | |
331 | Float_t ezm =0.; | |
332 | Float_t lim1 = hc->GetBinLowEdge(binmin); | |
333 | Float_t lim2 = hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax); | |
334 | Float_t widthSR=lim2-lim1; | |
335 | ||
336 | if(nPeaks ==1 && (lim2-lim1)<fWindowWidth){ | |
337 | Float_t c=(lim1+lim2)/2.; | |
338 | lim1=c-fWindowWidth/2.; | |
339 | lim2=c+fWindowWidth/2.; | |
340 | } | |
341 | Int_t niter = 0, ncontr=0; | |
342 | do { | |
343 | // symmetrization | |
344 | if(zm !=0.){ | |
345 | Float_t semilarg=TMath::Min((lim2-zm),(zm-lim1)); | |
346 | lim1=zm - semilarg; | |
347 | lim2=zm + semilarg; | |
348 | } | |
349 | ||
350 | zm=0.; | |
351 | ezm=0.; | |
352 | ncontr=0; | |
353 | for(Int_t i =0; i<points.GetEntriesFast(); i++){ | |
354 | AliITSZPoint* p=(AliITSZPoint*)points.UncheckedAt(i); | |
355 | if(p->GetZ()>lim1 && p->GetZ()<lim2){ | |
356 | Float_t deno = p->GetErrZ(); | |
357 | zm+=p->GetZ()/deno; | |
358 | ezm+=1./deno; | |
359 | ncontr++; | |
360 | } | |
361 | } | |
362 | if(ezm>0) { | |
363 | zm/=ezm; | |
364 | ezm=TMath::Sqrt(1./ezm); | |
365 | } | |
366 | niter++; | |
367 | } while(niter<10 && TMath::Abs((zm-lim1)-(lim2-zm))>fTolerance); | |
368 | if(nPeaks==2) ezm=widthSR; | |
369 | Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],zm}; | |
370 | Double_t covmatrix[6]={GetNominalCov()[0],0.,GetNominalCov()[2],0.,0.,ezm}; | |
371 | fCurrentVertex = new AliESDVertex(position,covmatrix,99999.,ncontr); | |
372 | fCurrentVertex->SetTitle("vertexer: Z"); | |
373 | fCurrentVertex->SetDispersion(fDiffPhiMax); | |
374 | fNoVertices=1; | |
375 | points.Clear(); | |
376 | if(fSearchForPileup && ncontr>fMinTrackletsForPilup){ | |
377 | Float_t secPeakPos; | |
378 | Int_t ncontr2=FindSecondPeak(fZCombc,binmin,binmax,secPeakPos); | |
379 | if(ncontr2>=fMinTrackletsForPilup){ | |
380 | fIsPileup=kTRUE; | |
381 | fNoVertices=2; | |
382 | fZpuv=secPeakPos; | |
383 | fNTrpuv=ncontr2; | |
384 | AliESDVertex secondVert(secPeakPos,0.1,ncontr2); | |
385 | fVertArray = new AliESDVertex[2]; | |
386 | fVertArray[0]=(*fCurrentVertex); | |
387 | fVertArray[1]=secondVert; | |
388 | } | |
389 | } | |
390 | if(fNoVertices==1){ | |
391 | fVertArray = new AliESDVertex[1]; | |
392 | fVertArray[0]=(*fCurrentVertex); | |
393 | } | |
394 | ||
395 | ResetHistograms(); | |
396 | return; | |
397 | } | |
398 | ||
399 | //_____________________________________________________________________ | |
400 | Int_t AliITSVertexerZD::FindSecondPeak(TH1F* h, Int_t binmin,Int_t binmax, Float_t& secPeakPos){ | |
401 | // Resets bin contents between binmin and binmax and then search | |
402 | // for a second peak position | |
403 | for(Int_t i=binmin-1;i<=binmax+1;i++){ | |
404 | h->SetBinContent(i,0.); | |
405 | } | |
406 | Int_t secPeakBin=h->GetMaximumBin(); | |
407 | secPeakPos=h->GetBinCenter(secPeakBin); | |
408 | Int_t secPeakCont=(Int_t)h->GetBinContent(secPeakBin); | |
409 | secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-1); | |
410 | secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+1); | |
411 | secPeakCont+=(Int_t)h->GetBinContent(secPeakBin-2); | |
412 | secPeakCont+=(Int_t)h->GetBinContent(secPeakBin+2); | |
413 | return secPeakCont; | |
414 | } | |
415 | ||
416 | //_____________________________________________________________________ | |
417 | void AliITSVertexerZD::ResetHistograms(){ | |
418 | // delete TH1 data members | |
419 | if(fZCombc)delete fZCombc; | |
420 | fZCombc = 0; | |
421 | } | |
422 | ||
423 | //________________________________________________________ | |
424 | void AliITSVertexerZD::PrintStatus() const { | |
425 | // Print current status | |
426 | cout <<"=======================================================\n"; | |
427 | cout <<" First layer first and last modules: "<<fFirstL1<<", "; | |
428 | cout <<fLastL1<<endl; | |
429 | cout <<" Second layer first and last modules: "<<fFirstL2<<", "; | |
430 | cout <<fLastL2<<endl; | |
431 | cout <<" Max Phi difference: "<<fDiffPhiMax<<endl; | |
432 | cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl; | |
433 | cout <<"Bin sizes for coarse z histos "<<fStepCoarse<<endl; | |
434 | cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl; | |
435 | if(fZCombc){ | |
436 | cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl; | |
437 | } | |
438 | else{ | |
439 | cout<<"fZCombc does not exist\n"; | |
440 | } | |
441 | ||
442 | cout <<"=======================================================\n"; | |
443 | } | |
444 |