]>
Commit | Line | Data |
---|---|---|
70c95f95 | 1 | /************************************************************************** |
2 | * Copyright(c) 2006-2008, 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 | **************************************************************************/ | |
70c95f95 | 15 | #include <TH3F.h> |
16 | #include <TTree.h> | |
6a0d56b8 | 17 | #include "AliESDVertex.h" |
cd706e57 | 18 | #include "AliLog.h" |
6a0d56b8 | 19 | #include "AliStrLine.h" |
20 | #include "AliTracker.h" | |
70c95f95 | 21 | #include "AliRunLoader.h" |
22 | #include "AliITSLoader.h" | |
70c95f95 | 23 | #include "AliITSDetTypeRec.h" |
24 | #include "AliITSRecPoint.h" | |
32e63e47 | 25 | #include "AliITSgeomTGeo.h" |
6a0d56b8 | 26 | #include "AliVertexerTracks.h" |
27 | #include "AliITSVertexer3D.h" | |
70c95f95 | 28 | ///////////////////////////////////////////////////////////////// |
29 | // this class implements a method to determine | |
30 | // the 3 coordinates of the primary vertex | |
31 | // for p-p collisions | |
32 | // It can be used successfully with Pb-Pb collisions | |
33 | //////////////////////////////////////////////////////////////// | |
34 | ||
35 | ClassImp(AliITSVertexer3D) | |
36 | ||
6a0d56b8 | 37 | /* $Id$ */ |
38 | ||
70c95f95 | 39 | //______________________________________________________________________ |
40 | AliITSVertexer3D::AliITSVertexer3D():AliITSVertexer(), | |
41 | fLines(), | |
42 | fVert3D(), | |
43 | fCoarseDiffPhiCut(0.), | |
05d1294c | 44 | fCoarseMaxRCut(0.), |
70c95f95 | 45 | fMaxRCut(0.), |
46 | fZCutDiamond(0.), | |
05d1294c | 47 | fMaxZCut(0.), |
70c95f95 | 48 | fDCAcut(0.), |
6a0d56b8 | 49 | fDiffPhiMax(0.), |
50 | fMeanPSelTrk(0.), | |
51 | fMeanPtSelTrk(0.) | |
52 | { | |
70c95f95 | 53 | // Default constructor |
54 | SetCoarseDiffPhiCut(); | |
05d1294c | 55 | SetCoarseMaxRCut(); |
70c95f95 | 56 | SetMaxRCut(); |
57 | SetZCutDiamond(); | |
05d1294c | 58 | SetMaxZCut(); |
70c95f95 | 59 | SetDCAcut(); |
60 | SetDiffPhiMax(); | |
6a0d56b8 | 61 | SetMeanPSelTracks(); |
62 | SetMeanPtSelTracks(); | |
70c95f95 | 63 | } |
64 | ||
65 | //______________________________________________________________________ | |
66 | AliITSVertexer3D::AliITSVertexer3D(TString fn): AliITSVertexer(fn), | |
f5f6da22 | 67 | fLines("AliStrLine",1000), |
70c95f95 | 68 | fVert3D(), |
05d1294c | 69 | fCoarseDiffPhiCut(0.), |
70 | fCoarseMaxRCut(0.), | |
70c95f95 | 71 | fMaxRCut(0.), |
72 | fZCutDiamond(0.), | |
05d1294c | 73 | fMaxZCut(0.), |
70c95f95 | 74 | fDCAcut(0.), |
6a0d56b8 | 75 | fDiffPhiMax(0.), |
76 | fMeanPSelTrk(0.), | |
77 | fMeanPtSelTrk(0.) | |
27167524 | 78 | { |
70c95f95 | 79 | // Standard constructor |
70c95f95 | 80 | SetCoarseDiffPhiCut(); |
05d1294c | 81 | SetCoarseMaxRCut(); |
70c95f95 | 82 | SetMaxRCut(); |
83 | SetZCutDiamond(); | |
05d1294c | 84 | SetMaxZCut(); |
70c95f95 | 85 | SetDCAcut(); |
86 | SetDiffPhiMax(); | |
6a0d56b8 | 87 | SetMeanPSelTracks(); |
88 | SetMeanPtSelTracks(); | |
70c95f95 | 89 | } |
90 | ||
91 | //______________________________________________________________________ | |
92 | AliITSVertexer3D::~AliITSVertexer3D() { | |
93 | // Destructor | |
f5f6da22 | 94 | fLines.Clear("C"); |
70c95f95 | 95 | } |
96 | ||
27167524 | 97 | //______________________________________________________________________ |
98 | void AliITSVertexer3D::ResetVert3D(){ | |
99 | // | |
100 | fVert3D.SetXv(0.); | |
101 | fVert3D.SetYv(0.); | |
102 | fVert3D.SetZv(0.); | |
103 | fVert3D.SetDispersion(0.); | |
104 | fVert3D.SetNContributors(0); | |
105 | } | |
70c95f95 | 106 | //______________________________________________________________________ |
107 | AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(Int_t evnumber){ | |
108 | // Defines the AliESDVertex for the current event | |
27167524 | 109 | ResetVert3D(); |
cd706e57 | 110 | AliDebug(1,Form("FindVertexForCurrentEvent - 3D - PROCESSING EVENT %d",evnumber)); |
f5f6da22 | 111 | fLines.Clear(); |
70c95f95 | 112 | |
113 | Int_t nolines = FindTracklets(evnumber,0); | |
114 | fCurrentVertex = 0; | |
115 | if(nolines<2)return fCurrentVertex; | |
05d1294c | 116 | Int_t rc=Prepare3DVertex(0); |
f5f6da22 | 117 | if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); |
27167524 | 118 | /* uncomment to debug |
119 | printf("Vertex found in first iteration:\n"); | |
120 | fVert3D.Print(); | |
121 | printf("Start second iteration\n"); | |
122 | end of debug lines */ | |
123 | if(fVert3D.GetNContributors()>0){ | |
f5f6da22 | 124 | fLines.Clear("C"); |
70c95f95 | 125 | nolines = FindTracklets(evnumber,1); |
d35ad08f | 126 | if(nolines>=2){ |
05d1294c | 127 | rc=Prepare3DVertex(1); |
f5f6da22 | 128 | if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); |
d35ad08f | 129 | } |
70c95f95 | 130 | } |
27167524 | 131 | /* uncomment to debug |
132 | printf("Vertex found in second iteration:\n"); | |
133 | fVert3D.Print(); | |
134 | end of debug lines */ | |
135 | ||
136 | Float_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv()); | |
137 | if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){ | |
6a0d56b8 | 138 | Double_t position[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()}; |
139 | Double_t covmatrix[6]; | |
140 | fVert3D.GetCovMatrix(covmatrix); | |
141 | Double_t chi2=99999.; | |
142 | Int_t nContr=fVert3D.GetNContributors(); | |
143 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
70c95f95 | 144 | fCurrentVertex->SetTitle("vertexer: 3D"); |
145 | fCurrentVertex->SetName("Vertex"); | |
27167524 | 146 | fCurrentVertex->SetDispersion(fVert3D.GetDispersion()); |
70c95f95 | 147 | } |
148 | FindMultiplicity(evnumber); | |
149 | return fCurrentVertex; | |
150 | } | |
151 | ||
152 | //______________________________________________________________________ | |
153 | Int_t AliITSVertexer3D::FindTracklets(Int_t evnumber, Int_t optCuts){ | |
154 | // All the possible combinations between recpoints on layer 1and 2 are | |
155 | // considered. Straight lines (=tracklets)are formed. | |
156 | // The tracklets are processed in Prepare3DVertex | |
157 | AliRunLoader *rl =AliRunLoader::GetRunLoader(); | |
158 | AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader"); | |
32e63e47 | 159 | // AliITSgeom* geom = itsLoader->GetITSgeom(); |
70c95f95 | 160 | itsLoader->LoadRecPoints(); |
161 | rl->GetEvent(evnumber); | |
162 | ||
163 | AliITSDetTypeRec detTypeRec; | |
164 | ||
165 | TTree *tR = itsLoader->TreeR(); | |
166 | detTypeRec.SetTreeAddressR(tR); | |
167 | TClonesArray *itsRec = 0; | |
6a0d56b8 | 168 | // lc1 and gc1 are local and global coordinates for layer 1 |
169 | // Float_t lc1[3]={0.,0.,0.}; | |
170 | Float_t gc1[3]={0.,0.,0.}; | |
70c95f95 | 171 | // lc2 and gc2 are local and global coordinates for layer 2 |
6a0d56b8 | 172 | // Float_t lc2[3]={0.,0.,0.}; |
173 | Float_t gc2[3]={0.,0.,0.}; | |
70c95f95 | 174 | |
175 | itsRec = detTypeRec.RecPoints(); | |
176 | TBranch *branch; | |
177 | branch = tR->GetBranch("ITSRecPoints"); | |
178 | ||
179 | // Set values for cuts | |
180 | Float_t xbeam=0., ybeam=0.; | |
05d1294c | 181 | Float_t zvert=0.; |
70c95f95 | 182 | Float_t deltaPhi=fCoarseDiffPhiCut; |
05d1294c | 183 | Float_t deltaR=fCoarseMaxRCut; |
184 | Float_t dZmax=fZCutDiamond; | |
70c95f95 | 185 | if(optCuts){ |
27167524 | 186 | xbeam=fVert3D.GetXv(); |
187 | ybeam=fVert3D.GetYv(); | |
188 | zvert=fVert3D.GetZv(); | |
70c95f95 | 189 | deltaPhi = fDiffPhiMax; |
05d1294c | 190 | deltaR=fMaxRCut; |
191 | dZmax=fMaxZCut; | |
70c95f95 | 192 | } |
70c95f95 | 193 | Int_t nrpL1 = 0; // number of rec points on layer 1 |
194 | Int_t nrpL2 = 0; // number of rec points on layer 2 | |
195 | ||
196 | // By default irstL1=0 and lastL1=79 | |
6a0d56b8 | 197 | Int_t firstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1); |
32e63e47 | 198 | Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1; |
6a0d56b8 | 199 | for(Int_t module= firstL1; module<=lastL1;module++){ // count number of recopints on layer 1 |
70c95f95 | 200 | branch->GetEvent(module); |
201 | nrpL1+= itsRec->GetEntries(); | |
202 | detTypeRec.ResetRecPoints(); | |
203 | } | |
6a0d56b8 | 204 | //By default firstL2=80 and lastL2=239 |
205 | Int_t firstL2 = AliITSgeomTGeo::GetModuleIndex(2,1,1); | |
32e63e47 | 206 | Int_t lastL2 = AliITSgeomTGeo::GetModuleIndex(3,1,1)-1; |
6a0d56b8 | 207 | for(Int_t module= firstL2; module<=lastL2;module++){ // count number of recopints on layer 2 |
70c95f95 | 208 | branch->GetEvent(module); |
209 | nrpL2+= itsRec->GetEntries(); | |
210 | detTypeRec.ResetRecPoints(); | |
211 | } | |
212 | if(nrpL1 == 0 || nrpL2 == 0){ | |
213 | itsLoader->UnloadRecPoints(); | |
214 | return -1; | |
215 | } | |
cd706e57 | 216 | AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2)); |
70c95f95 | 217 | |
218 | Double_t a[3]={xbeam,ybeam,0.}; | |
219 | Double_t b[3]={xbeam,ybeam,10.}; | |
220 | AliStrLine zeta(a,b,kTRUE); | |
6a0d56b8 | 221 | Float_t bField=AliTracker::GetBz()/10.; //T |
222 | SetMeanPPtSelTracks(bField); | |
70c95f95 | 223 | |
224 | Int_t nolines = 0; | |
225 | // Loop on modules of layer 1 | |
6a0d56b8 | 226 | for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1 |
27167524 | 227 | UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1 |
70c95f95 | 228 | branch->GetEvent(modul1); |
229 | Int_t nrecp1 = itsRec->GetEntries(); | |
f5f6da22 | 230 | static TClonesArray prpl1("AliITSRecPoint",nrecp1); |
231 | prpl1.SetOwner(); | |
70c95f95 | 232 | for(Int_t j=0;j<nrecp1;j++){ |
233 | AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j); | |
f5f6da22 | 234 | new(prpl1[j])AliITSRecPoint(*recp); |
70c95f95 | 235 | } |
236 | detTypeRec.ResetRecPoints(); | |
237 | for(Int_t j=0;j<nrecp1;j++){ | |
f5f6da22 | 238 | AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1.At(j); |
70c95f95 | 239 | // Local coordinates of this recpoint |
32e63e47 | 240 | /* |
6a0d56b8 | 241 | lc[0]=recp1->GetDetLocalX(); |
242 | lc[2]=recp1->GetDetLocalZ(); | |
32e63e47 | 243 | */ |
6a0d56b8 | 244 | recp1->GetGlobalXYZ(gc1); |
245 | Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam); | |
70c95f95 | 246 | if(phi1<0)phi1=2*TMath::Pi()+phi1; |
27167524 | 247 | for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){ |
248 | for(Int_t k=0;k<4;k++){ | |
249 | Int_t ladmod=fLadders[ladder-1]+ladl2; | |
32e63e47 | 250 | if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2); |
251 | Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1); | |
27167524 | 252 | branch->GetEvent(modul2); |
253 | Int_t nrecp2 = itsRec->GetEntries(); | |
254 | for(Int_t j2=0;j2<nrecp2;j2++){ | |
6a0d56b8 | 255 | AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2); |
32e63e47 | 256 | /* |
6a0d56b8 | 257 | lc2[0]=recp2->GetDetLocalX(); |
258 | lc2[2]=recp2->GetDetLocalZ(); | |
32e63e47 | 259 | */ |
6a0d56b8 | 260 | recp2->GetGlobalXYZ(gc2); |
27167524 | 261 | Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam); |
262 | if(phi2<0)phi2=2*TMath::Pi()+phi2; | |
263 | Double_t diff = TMath::Abs(phi2-phi1); | |
264 | if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; | |
265 | if(diff>deltaPhi)continue; | |
6a0d56b8 | 266 | AliStrLine line(gc1,gc2,kTRUE); |
27167524 | 267 | Double_t cp[3]; |
268 | Int_t retcode = line.Cross(&zeta,cp); | |
269 | if(retcode<0)continue; | |
270 | Double_t dca = line.GetDCA(&zeta); | |
271 | if(dca<0.) continue; | |
272 | if(dca>deltaR)continue; | |
273 | Double_t deltaZ=cp[2]-zvert; | |
274 | if(TMath::Abs(deltaZ)>dZmax)continue; | |
6a0d56b8 | 275 | |
6a0d56b8 | 276 | if(nolines == 0){ |
f5f6da22 | 277 | if(fLines.GetEntriesFast()>0)fLines.Clear("C"); |
6a0d56b8 | 278 | } |
279 | Float_t cov[6]; | |
280 | recp2->GetGlobalCov(cov); | |
281 | ||
282 | ||
283 | Float_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); | |
284 | Float_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); | |
285 | Float_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction | |
286 | ||
287 | Float_t curvErr=0; | |
288 | if(bField>0.00001){ | |
289 | Float_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm | |
290 | Float_t dRad=TMath::Sqrt(TMath::Power((gc1[0]-gc2[0]),2)+TMath::Power((gc1[1]-gc2[1]),2)); | |
291 | Float_t aux=dRad/2.+rad1; | |
292 | curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm | |
293 | } | |
294 | ||
295 | Float_t sigmasq[3]; | |
296 | sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor; | |
297 | sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor; | |
298 | sigmasq[2]=cov[5]*factor*factor; | |
299 | ||
300 | // Multiple scattering | |
301 | Float_t beta=1.; | |
302 | Float_t beta2=beta*beta; | |
303 | Float_t p2=fMeanPSelTrk*fMeanPSelTrk; | |
304 | Float_t rBP=GetPipeRadius(); | |
305 | Float_t dBP=0.08/35.3; // 800 um of Be | |
306 | Float_t dL1=0.01; //approx. 1% of radiation length | |
307 | Float_t theta2BP=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dBP); | |
308 | Float_t theta2L1=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dL1); | |
309 | Float_t thetaBP=TMath::Sqrt(theta2BP); | |
310 | Float_t thetaL1=TMath::Sqrt(theta2L1); | |
311 | // Float_t geomfac[3]; | |
312 | // geomfac[0]=sin(phi1)*sin(phi1); | |
313 | // geomfac[1]=cos(phi1)*cos(phi1); | |
314 | // Float_t tgth=(gc2[2]-gc1[2])/(rad2-rad1); | |
315 | // geomfac[2]=1+tgth*tgth; | |
316 | for(Int_t ico=0; ico<3;ico++){ | |
317 | // printf("Error on coord. %d due to cov matrix+curvErr=%f\n",ico,sigmasq[ico]); | |
318 | // // sigmasq[ico]+=rad1*rad1*geomfac[ico]*theta2L1/2; // multiple scattering in layer 1 | |
319 | // // sigmasq[ico]+=rBP*rBP*geomfac[ico]*theta2BP/2; // multiple scattering in beam pipe | |
320 | sigmasq[ico]+=TMath::Power(rad1*TMath::Tan(thetaL1),2)/3.; | |
321 | sigmasq[ico]+=TMath::Power(rBP*TMath::Tan(thetaBP),2)/3.; | |
322 | ||
323 | // printf("Multipl. scatt. contr %d = %f (LAY1), %f (BP)\n",ico,rad1*rad1*geomfac[ico]*theta2L1/2,rBP*rBP*geomfac[ico]*theta2BP/2); | |
324 | // printf("Total error on coord %d = %f\n",ico,sigmasq[ico]); | |
325 | } | |
326 | Float_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.}; | |
327 | if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0]; | |
328 | if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1]; | |
329 | if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2]; | |
f5f6da22 | 330 | new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE); |
6a0d56b8 | 331 | |
27167524 | 332 | } |
333 | detTypeRec.ResetRecPoints(); | |
70c95f95 | 334 | } |
70c95f95 | 335 | } |
336 | } | |
f5f6da22 | 337 | prpl1.Clear(); |
70c95f95 | 338 | } |
339 | if(nolines == 0)return -2; | |
340 | return nolines; | |
341 | } | |
342 | ||
343 | //______________________________________________________________________ | |
344 | void AliITSVertexer3D::FindVertices(){ | |
345 | // computes the vertices of the events in the range FirstEvent - LastEvent | |
346 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); | |
347 | AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader"); | |
348 | itsLoader->ReloadRecPoints(); | |
349 | for(Int_t i=fFirstEvent;i<=fLastEvent;i++){ | |
350 | rl->GetEvent(i); | |
351 | FindVertexForCurrentEvent(i); | |
352 | if(fCurrentVertex){ | |
353 | WriteCurrentVertex(); | |
354 | } | |
355 | } | |
356 | } | |
357 | ||
358 | ||
70c95f95 | 359 | //______________________________________________________________________ |
05d1294c | 360 | Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ |
70c95f95 | 361 | // Finds the 3D vertex information using tracklets |
362 | Int_t retcode = -1; | |
363 | ||
05d1294c | 364 | Float_t xbeam=0.; |
365 | Float_t ybeam=0.; | |
366 | Float_t zvert=0.; | |
367 | Float_t deltaR=fCoarseMaxRCut; | |
368 | Float_t dZmax=fZCutDiamond; | |
369 | if(optCuts){ | |
27167524 | 370 | xbeam=fVert3D.GetXv(); |
371 | ybeam=fVert3D.GetYv(); | |
372 | zvert=fVert3D.GetZv(); | |
05d1294c | 373 | deltaR=fMaxRCut; |
374 | dZmax=fMaxZCut; | |
375 | } | |
376 | ||
70c95f95 | 377 | Int_t nbr=50; |
05d1294c | 378 | Float_t rl=-fCoarseMaxRCut; |
379 | Float_t rh=fCoarseMaxRCut; | |
70c95f95 | 380 | Int_t nbz=100; |
381 | Float_t zl=-fZCutDiamond; | |
382 | Float_t zh=fZCutDiamond; | |
383 | Float_t binsizer=(rh-rl)/nbr; | |
384 | Float_t binsizez=(zh-zl)/nbz; | |
385 | TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh); | |
386 | ||
70c95f95 | 387 | // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed) |
f5f6da22 | 388 | Int_t *validate = new Int_t [fLines.GetEntriesFast()]; |
389 | for(Int_t i=0; i<fLines.GetEntriesFast();i++)validate[i]=0; | |
390 | for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){ | |
05d1294c | 391 | if(validate[i]==1)continue; |
f5f6da22 | 392 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); |
393 | for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){ | |
394 | AliStrLine *l2 = (AliStrLine*)fLines.At(j); | |
05d1294c | 395 | Float_t dca=l1->GetDCA(l2); |
396 | if(dca > fDCAcut || dca<0.00001) continue; | |
70c95f95 | 397 | Double_t point[3]; |
398 | Int_t retc = l1->Cross(l2,point); | |
399 | if(retc<0)continue; | |
6a0d56b8 | 400 | Double_t deltaZ=point[2]-zvert; |
401 | if(TMath::Abs(deltaZ)>dZmax)continue; | |
70c95f95 | 402 | Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]); |
05d1294c | 403 | if(rad>fCoarseMaxRCut)continue; |
404 | Double_t deltaX=point[0]-xbeam; | |
405 | Double_t deltaY=point[1]-ybeam; | |
05d1294c | 406 | Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY); |
05d1294c | 407 | if(raddist>deltaR)continue; |
27167524 | 408 | validate[i]=1; |
409 | validate[j]=1; | |
70c95f95 | 410 | h3d->Fill(point[0],point[1],point[2]); |
411 | } | |
412 | } | |
413 | ||
414 | ||
415 | ||
70c95f95 | 416 | Int_t numbtracklets=0; |
f5f6da22 | 417 | for(Int_t i=0; i<fLines.GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++; |
70c95f95 | 418 | if(numbtracklets<2){delete [] validate; delete h3d; return retcode; } |
419 | ||
f5f6da22 | 420 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
421 | if(validate[i]<1)fLines.RemoveAt(i); | |
70c95f95 | 422 | } |
f5f6da22 | 423 | fLines.Compress(); |
424 | AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast())); | |
70c95f95 | 425 | delete [] validate; |
426 | ||
05d1294c | 427 | |
428 | // finds peak in histo | |
70c95f95 | 429 | TAxis *xax = h3d->GetXaxis(); |
430 | TAxis *yax = h3d->GetYaxis(); | |
431 | TAxis *zax = h3d->GetZaxis(); | |
70c95f95 | 432 | Double_t peak[3]={0.,0.,0.}; |
433 | Float_t contref = 0.; | |
434 | for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){ | |
435 | Float_t xval = xax->GetBinCenter(i); | |
436 | for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){ | |
437 | Float_t yval = yax->GetBinCenter(j); | |
438 | for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){ | |
439 | Float_t bc = h3d->GetBinContent(i,j,k); | |
440 | Float_t zval = zax->GetBinCenter(k); | |
441 | if(bc>contref){ | |
442 | contref = bc; | |
443 | peak[2] = zval; | |
444 | peak[1] = yval; | |
445 | peak[0] = xval; | |
446 | } | |
447 | } | |
448 | } | |
449 | } | |
450 | delete h3d; | |
451 | ||
452 | // Second selection loop | |
453 | Float_t bs=(binsizer+binsizez)/2.; | |
f5f6da22 | 454 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
455 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); | |
456 | if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i); | |
70c95f95 | 457 | } |
f5f6da22 | 458 | fLines.Compress(); |
459 | AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast())); | |
70c95f95 | 460 | |
f5f6da22 | 461 | if(fLines.GetEntriesFast()>1){ |
27167524 | 462 | // find a first candidate for the primary vertex |
f5f6da22 | 463 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); |
70c95f95 | 464 | // make a further selection on tracklets based on this first candidate |
27167524 | 465 | fVert3D.GetXYZ(peak); |
466 | AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2])); | |
f5f6da22 | 467 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
468 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); | |
469 | if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i); | |
70c95f95 | 470 | } |
f5f6da22 | 471 | fLines.Compress(); |
472 | AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast())); | |
473 | if(fLines.GetEntriesFast()>1) retcode=0; // this new tracklet selection is used | |
27167524 | 474 | else retcode =1; // the previous tracklet selection will be used |
70c95f95 | 475 | } |
476 | else { | |
477 | retcode = 0; | |
478 | } | |
479 | return retcode; | |
480 | } | |
481 | ||
6a0d56b8 | 482 | //________________________________________________________ |
483 | void AliITSVertexer3D::SetMeanPPtSelTracks(Float_t fieldTesla){ | |
484 | // Sets mean values of P and Pt based on the field | |
485 | if(TMath::Abs(fieldTesla-0.5)<0.01){ | |
486 | SetMeanPSelTracks(0.885); | |
487 | SetMeanPtSelTracks(0.630); | |
488 | }else if(TMath::Abs(fieldTesla-0.4)<0.01){ | |
489 | SetMeanPSelTracks(0.805); | |
490 | SetMeanPtSelTracks(0.580); | |
491 | }else if(TMath::Abs(fieldTesla-0.2)<0.01){ | |
492 | SetMeanPSelTracks(0.740); | |
493 | SetMeanPtSelTracks(0.530); | |
494 | }else if(fieldTesla<0.00001){ | |
495 | SetMeanPSelTracks(0.730); | |
496 | SetMeanPtSelTracks(0.510); | |
497 | }else{ | |
498 | SetMeanPSelTracks(); | |
499 | SetMeanPtSelTracks(); | |
70c95f95 | 500 | } |
70c95f95 | 501 | } |
502 | ||
70c95f95 | 503 | |
504 | //________________________________________________________ | |
505 | void AliITSVertexer3D::PrintStatus() const { | |
506 | // Print current status | |
6a0d56b8 | 507 | printf("=======================================================\n"); |
508 | printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut); | |
509 | printf("Cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut); | |
510 | printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut); | |
511 | printf("Cut on diamond (Z) %f\n",fZCutDiamond); | |
512 | printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut); | |
513 | printf(" Max Phi difference: %f\n",fDiffPhiMax); | |
514 | printf("=======================================================\n"); | |
70c95f95 | 515 | |
70c95f95 | 516 | } |