]>
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 <TTree.h> |
6a0d56b8 | 16 | #include "AliESDVertex.h" |
cd706e57 | 17 | #include "AliLog.h" |
6a0d56b8 | 18 | #include "AliStrLine.h" |
19 | #include "AliTracker.h" | |
70c95f95 | 20 | #include "AliITSDetTypeRec.h" |
21 | #include "AliITSRecPoint.h" | |
32e63e47 | 22 | #include "AliITSgeomTGeo.h" |
6a0d56b8 | 23 | #include "AliVertexerTracks.h" |
24 | #include "AliITSVertexer3D.h" | |
0b7bac55 | 25 | #include "AliITSVertexerZ.h" |
70c95f95 | 26 | ///////////////////////////////////////////////////////////////// |
27 | // this class implements a method to determine | |
28 | // the 3 coordinates of the primary vertex | |
29 | // for p-p collisions | |
30 | // It can be used successfully with Pb-Pb collisions | |
31 | //////////////////////////////////////////////////////////////// | |
32 | ||
33 | ClassImp(AliITSVertexer3D) | |
34 | ||
6a0d56b8 | 35 | /* $Id$ */ |
36 | ||
70c95f95 | 37 | //______________________________________________________________________ |
38 | AliITSVertexer3D::AliITSVertexer3D():AliITSVertexer(), | |
308c2f7c | 39 | fLines("AliStrLine",1000), |
70c95f95 | 40 | fVert3D(), |
41 | fCoarseDiffPhiCut(0.), | |
05d1294c | 42 | fCoarseMaxRCut(0.), |
70c95f95 | 43 | fMaxRCut(0.), |
44 | fZCutDiamond(0.), | |
05d1294c | 45 | fMaxZCut(0.), |
70c95f95 | 46 | fDCAcut(0.), |
6a0d56b8 | 47 | fDiffPhiMax(0.), |
48 | fMeanPSelTrk(0.), | |
cd1c2af1 | 49 | fMeanPtSelTrk(0.), |
50 | fUsedCluster(kMaxCluPerMod*kNSPDMod), | |
51 | fZHisto(0), | |
52 | fDCAforPileup(0.), | |
53 | fPileupAlgo(0) | |
6a0d56b8 | 54 | { |
70c95f95 | 55 | // Default constructor |
56 | SetCoarseDiffPhiCut(); | |
05d1294c | 57 | SetCoarseMaxRCut(); |
70c95f95 | 58 | SetMaxRCut(); |
59 | SetZCutDiamond(); | |
05d1294c | 60 | SetMaxZCut(); |
7203e11a | 61 | SetDCACut(); |
70c95f95 | 62 | SetDiffPhiMax(); |
6a0d56b8 | 63 | SetMeanPSelTracks(); |
64 | SetMeanPtSelTracks(); | |
cd1c2af1 | 65 | SetMinDCAforPileup(); |
66 | SetPileupAlgo(); | |
67 | Float_t binsize=0.02; // default 200 micron | |
68 | Int_t nbins=static_cast<Int_t>(1+2*fZCutDiamond/binsize); | |
69 | fZHisto=new TH1F("hz","",nbins,-fZCutDiamond,-fZCutDiamond+binsize*nbins); | |
70c95f95 | 70 | } |
71 | ||
70c95f95 | 72 | //______________________________________________________________________ |
73 | AliITSVertexer3D::~AliITSVertexer3D() { | |
74 | // Destructor | |
f5f6da22 | 75 | fLines.Clear("C"); |
cd1c2af1 | 76 | if(fZHisto) delete fZHisto; |
70c95f95 | 77 | } |
78 | ||
27167524 | 79 | //______________________________________________________________________ |
80 | void AliITSVertexer3D::ResetVert3D(){ | |
81 | // | |
82 | fVert3D.SetXv(0.); | |
83 | fVert3D.SetYv(0.); | |
84 | fVert3D.SetZv(0.); | |
85 | fVert3D.SetDispersion(0.); | |
86 | fVert3D.SetNContributors(0); | |
cd1c2af1 | 87 | fIsPileup=kFALSE; |
88 | fNTrpuv=-2; | |
89 | fZpuv=-99999.; | |
90 | fUsedCluster.ResetAllBits(0); | |
27167524 | 91 | } |
70c95f95 | 92 | //______________________________________________________________________ |
308c2f7c | 93 | AliESDVertex* AliITSVertexer3D::FindVertexForCurrentEvent(TTree *itsClusterTree){ |
70c95f95 | 94 | // Defines the AliESDVertex for the current event |
27167524 | 95 | ResetVert3D(); |
308c2f7c | 96 | AliDebug(1,"FindVertexForCurrentEvent - 3D - PROCESSING NEXT EVENT"); |
f5f6da22 | 97 | fLines.Clear(); |
70c95f95 | 98 | |
308c2f7c | 99 | Int_t nolines = FindTracklets(itsClusterTree,0); |
70c95f95 | 100 | fCurrentVertex = 0; |
0b7bac55 | 101 | if(nolines>=2){ |
102 | Int_t rc=Prepare3DVertex(0); | |
103 | if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); | |
27167524 | 104 | /* uncomment to debug |
0b7bac55 | 105 | printf("Vertex found in first iteration:\n"); |
106 | fVert3D.Print(); | |
107 | printf("Start second iteration\n"); | |
27167524 | 108 | end of debug lines */ |
0b7bac55 | 109 | if(fVert3D.GetNContributors()>0){ |
110 | fLines.Clear("C"); | |
111 | nolines = FindTracklets(itsClusterTree,1); | |
112 | if(nolines>=2){ | |
113 | rc=Prepare3DVertex(1); | |
114 | if(rc==0) fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); | |
115 | } | |
d35ad08f | 116 | } |
27167524 | 117 | /* uncomment to debug |
0b7bac55 | 118 | printf("Vertex found in second iteration:\n"); |
119 | fVert3D.Print(); | |
120 | end of debug lines */ | |
27167524 | 121 | |
0b7bac55 | 122 | Float_t vRadius=TMath::Sqrt(fVert3D.GetXv()*fVert3D.GetXv()+fVert3D.GetYv()*fVert3D.GetYv()); |
123 | if(vRadius<GetPipeRadius() && fVert3D.GetNContributors()>0){ | |
124 | Double_t position[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()}; | |
125 | Double_t covmatrix[6]; | |
126 | fVert3D.GetCovMatrix(covmatrix); | |
127 | Double_t chi2=99999.; | |
128 | Int_t nContr=fVert3D.GetNContributors(); | |
129 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
130 | fCurrentVertex->SetTitle("vertexer: 3D"); | |
131 | fCurrentVertex->SetName("SPDVertex3D"); | |
132 | fCurrentVertex->SetDispersion(fVert3D.GetDispersion()); | |
cd1c2af1 | 133 | |
134 | switch(fPileupAlgo){ | |
135 | case 0: PileupFromZ(); break; | |
136 | case 1: FindOther3DVertices(itsClusterTree); break; | |
137 | default: AliError("Wrong pileup algorithm"); break; | |
138 | } | |
0b7bac55 | 139 | } |
140 | } | |
141 | if(!fCurrentVertex){ | |
142 | AliITSVertexerZ vertz(GetNominalPos()[0],GetNominalPos()[1]); | |
b96ee725 | 143 | vertz.SetDetTypeRec(GetDetTypeRec()); |
0b7bac55 | 144 | AliDebug(1,"Call Vertexer Z\n"); |
cc282e54 | 145 | vertz.SetLowLimit(-fZCutDiamond); |
146 | vertz.SetHighLimit(fZCutDiamond); | |
0b7bac55 | 147 | AliESDVertex* vtxz = vertz.FindVertexForCurrentEvent(itsClusterTree); |
148 | if(vtxz){ | |
149 | Double_t position[3]={GetNominalPos()[0],GetNominalPos()[1],vtxz->GetZv()}; | |
150 | Double_t covmatrix[6]; | |
151 | vtxz->GetCovMatrix(covmatrix); | |
152 | Double_t chi2=99999.; | |
153 | Int_t nContr=vtxz->GetNContributors(); | |
154 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nContr); | |
155 | fCurrentVertex->SetTitle("vertexer: Z"); | |
156 | fCurrentVertex->SetName("SPDVertexZ"); | |
157 | delete vtxz; | |
158 | // printf("Vertexer Z success\n"); | |
159 | } | |
160 | ||
70c95f95 | 161 | } |
308c2f7c | 162 | FindMultiplicity(itsClusterTree); |
70c95f95 | 163 | return fCurrentVertex; |
164 | } | |
165 | ||
166 | //______________________________________________________________________ | |
308c2f7c | 167 | Int_t AliITSVertexer3D::FindTracklets(TTree *itsClusterTree, Int_t optCuts){ |
70c95f95 | 168 | // All the possible combinations between recpoints on layer 1and 2 are |
169 | // considered. Straight lines (=tracklets)are formed. | |
170 | // The tracklets are processed in Prepare3DVertex | |
70c95f95 | 171 | |
308c2f7c | 172 | TTree *tR = itsClusterTree; |
06a7cbee | 173 | fDetTypeRec->SetTreeAddressR(tR); |
70c95f95 | 174 | TClonesArray *itsRec = 0; |
cd1c2af1 | 175 | if(optCuts==0) fZHisto->Reset(); |
6a0d56b8 | 176 | // lc1 and gc1 are local and global coordinates for layer 1 |
177 | // Float_t lc1[3]={0.,0.,0.}; | |
178 | Float_t gc1[3]={0.,0.,0.}; | |
70c95f95 | 179 | // lc2 and gc2 are local and global coordinates for layer 2 |
6a0d56b8 | 180 | // Float_t lc2[3]={0.,0.,0.}; |
181 | Float_t gc2[3]={0.,0.,0.}; | |
70c95f95 | 182 | |
06a7cbee | 183 | itsRec = fDetTypeRec->RecPoints(); |
70c95f95 | 184 | TBranch *branch; |
185 | branch = tR->GetBranch("ITSRecPoints"); | |
186 | ||
187 | // Set values for cuts | |
ecce5370 | 188 | Float_t xbeam=GetNominalPos()[0]; |
189 | Float_t ybeam=GetNominalPos()[1]; | |
05d1294c | 190 | Float_t zvert=0.; |
70c95f95 | 191 | Float_t deltaPhi=fCoarseDiffPhiCut; |
05d1294c | 192 | Float_t deltaR=fCoarseMaxRCut; |
193 | Float_t dZmax=fZCutDiamond; | |
cd1c2af1 | 194 | if(optCuts==1){ |
27167524 | 195 | xbeam=fVert3D.GetXv(); |
196 | ybeam=fVert3D.GetYv(); | |
197 | zvert=fVert3D.GetZv(); | |
70c95f95 | 198 | deltaPhi = fDiffPhiMax; |
05d1294c | 199 | deltaR=fMaxRCut; |
200 | dZmax=fMaxZCut; | |
cd1c2af1 | 201 | }else if(optCuts==2){ |
202 | xbeam=fVert3D.GetXv(); | |
203 | ybeam=fVert3D.GetYv(); | |
204 | deltaPhi = fDiffPhiMax; | |
205 | deltaR=fMaxRCut; | |
70c95f95 | 206 | } |
cd1c2af1 | 207 | |
70c95f95 | 208 | Int_t nrpL1 = 0; // number of rec points on layer 1 |
209 | Int_t nrpL2 = 0; // number of rec points on layer 2 | |
210 | ||
211 | // By default irstL1=0 and lastL1=79 | |
6a0d56b8 | 212 | Int_t firstL1 = AliITSgeomTGeo::GetModuleIndex(1,1,1); |
32e63e47 | 213 | Int_t lastL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1)-1; |
6a0d56b8 | 214 | for(Int_t module= firstL1; module<=lastL1;module++){ // count number of recopints on layer 1 |
70c95f95 | 215 | branch->GetEvent(module); |
216 | nrpL1+= itsRec->GetEntries(); | |
06a7cbee | 217 | fDetTypeRec->ResetRecPoints(); |
70c95f95 | 218 | } |
6a0d56b8 | 219 | //By default firstL2=80 and lastL2=239 |
220 | Int_t firstL2 = AliITSgeomTGeo::GetModuleIndex(2,1,1); | |
32e63e47 | 221 | Int_t lastL2 = AliITSgeomTGeo::GetModuleIndex(3,1,1)-1; |
6a0d56b8 | 222 | for(Int_t module= firstL2; module<=lastL2;module++){ // count number of recopints on layer 2 |
70c95f95 | 223 | branch->GetEvent(module); |
224 | nrpL2+= itsRec->GetEntries(); | |
06a7cbee | 225 | fDetTypeRec->ResetRecPoints(); |
70c95f95 | 226 | } |
227 | if(nrpL1 == 0 || nrpL2 == 0){ | |
70c95f95 | 228 | return -1; |
229 | } | |
cd706e57 | 230 | AliDebug(1,Form("RecPoints on Layer 1,2 = %d, %d\n",nrpL1,nrpL2)); |
70c95f95 | 231 | |
232 | Double_t a[3]={xbeam,ybeam,0.}; | |
233 | Double_t b[3]={xbeam,ybeam,10.}; | |
234 | AliStrLine zeta(a,b,kTRUE); | |
6a0d56b8 | 235 | Float_t bField=AliTracker::GetBz()/10.; //T |
236 | SetMeanPPtSelTracks(bField); | |
70c95f95 | 237 | |
238 | Int_t nolines = 0; | |
239 | // Loop on modules of layer 1 | |
6a0d56b8 | 240 | for(Int_t modul1= firstL1; modul1<=lastL1;modul1++){ // Loop on modules of layer 1 |
8c42830a | 241 | if(!fUseModule[modul1]) continue; |
27167524 | 242 | UShort_t ladder=int(modul1/4)+1; // ladders are numbered starting from 1 |
70c95f95 | 243 | branch->GetEvent(modul1); |
244 | Int_t nrecp1 = itsRec->GetEntries(); | |
f5f6da22 | 245 | static TClonesArray prpl1("AliITSRecPoint",nrecp1); |
246 | prpl1.SetOwner(); | |
70c95f95 | 247 | for(Int_t j=0;j<nrecp1;j++){ |
248 | AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j); | |
f5f6da22 | 249 | new(prpl1[j])AliITSRecPoint(*recp); |
70c95f95 | 250 | } |
06a7cbee | 251 | fDetTypeRec->ResetRecPoints(); |
70c95f95 | 252 | for(Int_t j=0;j<nrecp1;j++){ |
cd1c2af1 | 253 | if(j>kMaxCluPerMod) continue; |
254 | UShort_t idClu1=modul1*kMaxCluPerMod+j; | |
255 | if(fUsedCluster.TestBitNumber(idClu1)) continue; | |
f5f6da22 | 256 | AliITSRecPoint *recp1 = (AliITSRecPoint*)prpl1.At(j); |
70c95f95 | 257 | // Local coordinates of this recpoint |
32e63e47 | 258 | /* |
6a0d56b8 | 259 | lc[0]=recp1->GetDetLocalX(); |
260 | lc[2]=recp1->GetDetLocalZ(); | |
32e63e47 | 261 | */ |
6a0d56b8 | 262 | recp1->GetGlobalXYZ(gc1); |
263 | Double_t phi1 = TMath::ATan2(gc1[1]-ybeam,gc1[0]-xbeam); | |
70c95f95 | 264 | if(phi1<0)phi1=2*TMath::Pi()+phi1; |
27167524 | 265 | for(Int_t ladl2=0 ; ladl2<fLadOnLay2*2+1;ladl2++){ |
266 | for(Int_t k=0;k<4;k++){ | |
267 | Int_t ladmod=fLadders[ladder-1]+ladl2; | |
32e63e47 | 268 | if(ladmod>AliITSgeomTGeo::GetNLadders(2)) ladmod=ladmod-AliITSgeomTGeo::GetNLadders(2); |
269 | Int_t modul2=AliITSgeomTGeo::GetModuleIndex(2,ladmod,k+1); | |
8c42830a | 270 | if(!fUseModule[modul2]) continue; |
27167524 | 271 | branch->GetEvent(modul2); |
272 | Int_t nrecp2 = itsRec->GetEntries(); | |
273 | for(Int_t j2=0;j2<nrecp2;j2++){ | |
cd1c2af1 | 274 | if(j2>kMaxCluPerMod) continue; |
275 | UShort_t idClu2=modul2*kMaxCluPerMod+j2; | |
276 | if(fUsedCluster.TestBitNumber(idClu2)) continue; | |
6a0d56b8 | 277 | AliITSRecPoint *recp2 = (AliITSRecPoint*)itsRec->At(j2); |
32e63e47 | 278 | /* |
6a0d56b8 | 279 | lc2[0]=recp2->GetDetLocalX(); |
280 | lc2[2]=recp2->GetDetLocalZ(); | |
32e63e47 | 281 | */ |
6a0d56b8 | 282 | recp2->GetGlobalXYZ(gc2); |
27167524 | 283 | Double_t phi2 = TMath::ATan2(gc2[1]-ybeam,gc2[0]-xbeam); |
284 | if(phi2<0)phi2=2*TMath::Pi()+phi2; | |
285 | Double_t diff = TMath::Abs(phi2-phi1); | |
286 | if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; | |
cd1c2af1 | 287 | if(optCuts==0 && diff<fDiffPhiMax){ |
288 | Float_t r1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); | |
289 | Float_t zc1=gc1[2]; | |
290 | Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); | |
291 | Float_t zc2=gc2[2]; | |
292 | Float_t zr0=(r2*zc1-r1*zc2)/(r2-r1); //Z @ null radius | |
293 | fZHisto->Fill(zr0); | |
294 | } | |
27167524 | 295 | if(diff>deltaPhi)continue; |
6a0d56b8 | 296 | AliStrLine line(gc1,gc2,kTRUE); |
27167524 | 297 | Double_t cp[3]; |
298 | Int_t retcode = line.Cross(&zeta,cp); | |
299 | if(retcode<0)continue; | |
300 | Double_t dca = line.GetDCA(&zeta); | |
301 | if(dca<0.) continue; | |
302 | if(dca>deltaR)continue; | |
303 | Double_t deltaZ=cp[2]-zvert; | |
304 | if(TMath::Abs(deltaZ)>dZmax)continue; | |
6a0d56b8 | 305 | |
6a0d56b8 | 306 | if(nolines == 0){ |
f5f6da22 | 307 | if(fLines.GetEntriesFast()>0)fLines.Clear("C"); |
6a0d56b8 | 308 | } |
309 | Float_t cov[6]; | |
310 | recp2->GetGlobalCov(cov); | |
311 | ||
cd1c2af1 | 312 | |
6a0d56b8 | 313 | Float_t rad1=TMath::Sqrt(gc1[0]*gc1[0]+gc1[1]*gc1[1]); |
314 | Float_t rad2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); | |
315 | Float_t factor=(rad1+rad2)/(rad2-rad1); //factor to account for error on tracklet direction | |
316 | ||
317 | Float_t curvErr=0; | |
318 | if(bField>0.00001){ | |
319 | Float_t curvRadius=fMeanPtSelTrk/(0.3*bField)*100; //cm | |
320 | Float_t dRad=TMath::Sqrt(TMath::Power((gc1[0]-gc2[0]),2)+TMath::Power((gc1[1]-gc2[1]),2)); | |
321 | Float_t aux=dRad/2.+rad1; | |
322 | curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm | |
323 | } | |
324 | ||
325 | Float_t sigmasq[3]; | |
326 | sigmasq[0]=(cov[0]+curvErr*curvErr/2.)*factor*factor; | |
327 | sigmasq[1]=(cov[3]+curvErr*curvErr/2.)*factor*factor; | |
328 | sigmasq[2]=cov[5]*factor*factor; | |
329 | ||
330 | // Multiple scattering | |
331 | Float_t beta=1.; | |
332 | Float_t beta2=beta*beta; | |
333 | Float_t p2=fMeanPSelTrk*fMeanPSelTrk; | |
334 | Float_t rBP=GetPipeRadius(); | |
335 | Float_t dBP=0.08/35.3; // 800 um of Be | |
336 | Float_t dL1=0.01; //approx. 1% of radiation length | |
337 | Float_t theta2BP=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dBP); | |
338 | Float_t theta2L1=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dL1); | |
339 | Float_t thetaBP=TMath::Sqrt(theta2BP); | |
340 | Float_t thetaL1=TMath::Sqrt(theta2L1); | |
341 | // Float_t geomfac[3]; | |
342 | // geomfac[0]=sin(phi1)*sin(phi1); | |
343 | // geomfac[1]=cos(phi1)*cos(phi1); | |
344 | // Float_t tgth=(gc2[2]-gc1[2])/(rad2-rad1); | |
345 | // geomfac[2]=1+tgth*tgth; | |
346 | for(Int_t ico=0; ico<3;ico++){ | |
347 | // printf("Error on coord. %d due to cov matrix+curvErr=%f\n",ico,sigmasq[ico]); | |
348 | // // sigmasq[ico]+=rad1*rad1*geomfac[ico]*theta2L1/2; // multiple scattering in layer 1 | |
349 | // // sigmasq[ico]+=rBP*rBP*geomfac[ico]*theta2BP/2; // multiple scattering in beam pipe | |
350 | sigmasq[ico]+=TMath::Power(rad1*TMath::Tan(thetaL1),2)/3.; | |
351 | sigmasq[ico]+=TMath::Power(rBP*TMath::Tan(thetaBP),2)/3.; | |
352 | ||
353 | // printf("Multipl. scatt. contr %d = %f (LAY1), %f (BP)\n",ico,rad1*rad1*geomfac[ico]*theta2L1/2,rBP*rBP*geomfac[ico]*theta2BP/2); | |
354 | // printf("Total error on coord %d = %f\n",ico,sigmasq[ico]); | |
355 | } | |
356 | Float_t wmat[9]={1.,0.,0.,0.,1.,0.,0.,0.,1.}; | |
357 | if(sigmasq[0]!=0.) wmat[0]=1./sigmasq[0]; | |
358 | if(sigmasq[1]!=0.) wmat[4]=1./sigmasq[1]; | |
359 | if(sigmasq[2]!=0.) wmat[8]=1./sigmasq[2]; | |
cd1c2af1 | 360 | new(fLines[nolines++])AliStrLine(gc1,sigmasq,wmat,gc2,kTRUE,idClu1,idClu2); |
6a0d56b8 | 361 | |
27167524 | 362 | } |
06a7cbee | 363 | fDetTypeRec->ResetRecPoints(); |
70c95f95 | 364 | } |
70c95f95 | 365 | } |
366 | } | |
f5f6da22 | 367 | prpl1.Clear(); |
70c95f95 | 368 | } |
369 | if(nolines == 0)return -2; | |
370 | return nolines; | |
371 | } | |
372 | ||
70c95f95 | 373 | //______________________________________________________________________ |
05d1294c | 374 | Int_t AliITSVertexer3D::Prepare3DVertex(Int_t optCuts){ |
70c95f95 | 375 | // Finds the 3D vertex information using tracklets |
376 | Int_t retcode = -1; | |
377 | ||
ecce5370 | 378 | Float_t xbeam=GetNominalPos()[0]; |
379 | Float_t ybeam=GetNominalPos()[1]; | |
05d1294c | 380 | Float_t zvert=0.; |
381 | Float_t deltaR=fCoarseMaxRCut; | |
382 | Float_t dZmax=fZCutDiamond; | |
cd1c2af1 | 383 | if(optCuts==1){ |
27167524 | 384 | xbeam=fVert3D.GetXv(); |
385 | ybeam=fVert3D.GetYv(); | |
386 | zvert=fVert3D.GetZv(); | |
05d1294c | 387 | deltaR=fMaxRCut; |
388 | dZmax=fMaxZCut; | |
cd1c2af1 | 389 | }else if(optCuts==2){ |
390 | xbeam=fVert3D.GetXv(); | |
391 | ybeam=fVert3D.GetYv(); | |
392 | deltaR=fMaxRCut; | |
05d1294c | 393 | } |
394 | ||
70c95f95 | 395 | Int_t nbr=50; |
05d1294c | 396 | Float_t rl=-fCoarseMaxRCut; |
397 | Float_t rh=fCoarseMaxRCut; | |
70c95f95 | 398 | Int_t nbz=100; |
399 | Float_t zl=-fZCutDiamond; | |
400 | Float_t zh=fZCutDiamond; | |
401 | Float_t binsizer=(rh-rl)/nbr; | |
402 | Float_t binsizez=(zh-zl)/nbz; | |
403 | TH3F *h3d = new TH3F("h3d","xyz distribution",nbr,rl,rh,nbr,rl,rh,nbz,zl,zh); | |
7340864d | 404 | Int_t nbrcs=25; |
405 | Int_t nbzcs=50; | |
406 | TH3F *h3dcs = new TH3F("h3dcs","xyz distribution",nbrcs,rl,rh,nbrcs,rl,rh,nbzcs,zl,zh); | |
70c95f95 | 407 | |
70c95f95 | 408 | // cleanup of the TCLonesArray of tracklets (i.e. fakes are removed) |
f5f6da22 | 409 | Int_t *validate = new Int_t [fLines.GetEntriesFast()]; |
410 | for(Int_t i=0; i<fLines.GetEntriesFast();i++)validate[i]=0; | |
411 | for(Int_t i=0; i<fLines.GetEntriesFast()-1;i++){ | |
f5f6da22 | 412 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); |
413 | for(Int_t j=i+1;j<fLines.GetEntriesFast();j++){ | |
414 | AliStrLine *l2 = (AliStrLine*)fLines.At(j); | |
05d1294c | 415 | Float_t dca=l1->GetDCA(l2); |
416 | if(dca > fDCAcut || dca<0.00001) continue; | |
70c95f95 | 417 | Double_t point[3]; |
418 | Int_t retc = l1->Cross(l2,point); | |
419 | if(retc<0)continue; | |
6a0d56b8 | 420 | Double_t deltaZ=point[2]-zvert; |
421 | if(TMath::Abs(deltaZ)>dZmax)continue; | |
70c95f95 | 422 | Double_t rad=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]); |
05d1294c | 423 | if(rad>fCoarseMaxRCut)continue; |
424 | Double_t deltaX=point[0]-xbeam; | |
425 | Double_t deltaY=point[1]-ybeam; | |
05d1294c | 426 | Double_t raddist=TMath::Sqrt(deltaX*deltaX+deltaY*deltaY); |
05d1294c | 427 | if(raddist>deltaR)continue; |
27167524 | 428 | validate[i]=1; |
429 | validate[j]=1; | |
70c95f95 | 430 | h3d->Fill(point[0],point[1],point[2]); |
7340864d | 431 | h3dcs->Fill(point[0],point[1],point[2]); |
70c95f95 | 432 | } |
433 | } | |
434 | ||
435 | ||
436 | ||
70c95f95 | 437 | Int_t numbtracklets=0; |
f5f6da22 | 438 | for(Int_t i=0; i<fLines.GetEntriesFast();i++)if(validate[i]>=1)numbtracklets++; |
e6fe443d | 439 | if(numbtracklets<2){delete [] validate; delete h3d; delete h3dcs; return retcode; } |
70c95f95 | 440 | |
f5f6da22 | 441 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
442 | if(validate[i]<1)fLines.RemoveAt(i); | |
70c95f95 | 443 | } |
f5f6da22 | 444 | fLines.Compress(); |
445 | AliDebug(1,Form("Number of tracklets (after compress)%d ",fLines.GetEntriesFast())); | |
70c95f95 | 446 | delete [] validate; |
447 | ||
7340864d | 448 | // Find peaks in histos |
05d1294c | 449 | |
70c95f95 | 450 | Double_t peak[3]={0.,0.,0.}; |
7340864d | 451 | Int_t ntrkl,ntimes; |
452 | FindPeaks(h3d,peak,ntrkl,ntimes); | |
70c95f95 | 453 | delete h3d; |
454 | ||
7340864d | 455 | if(optCuts==0 && ntrkl<=2){ |
456 | ntrkl=0; | |
457 | ntimes=0; | |
458 | FindPeaks(h3dcs,peak,ntrkl,ntimes); | |
459 | binsizer=(rh-rl)/nbrcs; | |
460 | binsizez=(zh-zl)/nbzcs; | |
461 | if(ntrkl==1 || ntimes>1){delete h3dcs; return retcode;} | |
462 | } | |
463 | delete h3dcs; | |
464 | ||
70c95f95 | 465 | // Second selection loop |
7340864d | 466 | |
70c95f95 | 467 | Float_t bs=(binsizer+binsizez)/2.; |
f5f6da22 | 468 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
469 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); | |
470 | if(l1->GetDistFromPoint(peak)>2.5*bs)fLines.RemoveAt(i); | |
70c95f95 | 471 | } |
f5f6da22 | 472 | fLines.Compress(); |
473 | AliDebug(1,Form("Number of tracklets (after 2nd compression) %d",fLines.GetEntriesFast())); | |
70c95f95 | 474 | |
f5f6da22 | 475 | if(fLines.GetEntriesFast()>1){ |
27167524 | 476 | // find a first candidate for the primary vertex |
f5f6da22 | 477 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); |
70c95f95 | 478 | // make a further selection on tracklets based on this first candidate |
27167524 | 479 | fVert3D.GetXYZ(peak); |
480 | AliDebug(1,Form("FIRST V candidate: %f ; %f ; %f",peak[0],peak[1],peak[2])); | |
f5f6da22 | 481 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ |
482 | AliStrLine *l1 = (AliStrLine*)fLines.At(i); | |
483 | if(l1->GetDistFromPoint(peak)> fDCAcut)fLines.RemoveAt(i); | |
70c95f95 | 484 | } |
f5f6da22 | 485 | fLines.Compress(); |
486 | AliDebug(1,Form("Number of tracklets (after 3rd compression) %d",fLines.GetEntriesFast())); | |
487 | if(fLines.GetEntriesFast()>1) retcode=0; // this new tracklet selection is used | |
27167524 | 488 | else retcode =1; // the previous tracklet selection will be used |
70c95f95 | 489 | } |
490 | else { | |
491 | retcode = 0; | |
492 | } | |
493 | return retcode; | |
494 | } | |
495 | ||
6a0d56b8 | 496 | //________________________________________________________ |
497 | void AliITSVertexer3D::SetMeanPPtSelTracks(Float_t fieldTesla){ | |
498 | // Sets mean values of P and Pt based on the field | |
499 | if(TMath::Abs(fieldTesla-0.5)<0.01){ | |
500 | SetMeanPSelTracks(0.885); | |
501 | SetMeanPtSelTracks(0.630); | |
502 | }else if(TMath::Abs(fieldTesla-0.4)<0.01){ | |
503 | SetMeanPSelTracks(0.805); | |
504 | SetMeanPtSelTracks(0.580); | |
505 | }else if(TMath::Abs(fieldTesla-0.2)<0.01){ | |
506 | SetMeanPSelTracks(0.740); | |
507 | SetMeanPtSelTracks(0.530); | |
508 | }else if(fieldTesla<0.00001){ | |
509 | SetMeanPSelTracks(0.730); | |
510 | SetMeanPtSelTracks(0.510); | |
511 | }else{ | |
512 | SetMeanPSelTracks(); | |
513 | SetMeanPtSelTracks(); | |
70c95f95 | 514 | } |
70c95f95 | 515 | } |
516 | ||
7340864d | 517 | //________________________________________________________ |
518 | void AliITSVertexer3D::FindPeaks(TH3F* histo, Double_t *peak, Int_t &nOfTracklets, Int_t &nOfTimes){ | |
519 | // Finds bin with max contents in 3D histo of tracket intersections | |
520 | TAxis *xax = histo->GetXaxis(); | |
521 | TAxis *yax = histo->GetYaxis(); | |
522 | TAxis *zax = histo->GetZaxis(); | |
523 | peak[0]=0.; | |
524 | peak[1]=0.; | |
525 | peak[2]=0.; | |
526 | nOfTracklets = 0; | |
527 | nOfTimes=0; | |
528 | for(Int_t i=xax->GetFirst();i<=xax->GetLast();i++){ | |
529 | Float_t xval = xax->GetBinCenter(i); | |
530 | for(Int_t j=yax->GetFirst();j<=yax->GetLast();j++){ | |
531 | Float_t yval = yax->GetBinCenter(j); | |
532 | for(Int_t k=zax->GetFirst();k<=zax->GetLast();k++){ | |
533 | Float_t zval = zax->GetBinCenter(k); | |
534 | Int_t bc =(Int_t)histo->GetBinContent(i,j,k); | |
535 | if(bc>nOfTracklets){ | |
536 | nOfTracklets = bc; | |
537 | peak[2] = zval; | |
538 | peak[1] = yval; | |
539 | peak[0] = xval; | |
540 | nOfTimes = 1; | |
06a7cbee | 541 | }else if(bc==nOfTracklets){ |
7340864d | 542 | nOfTimes++; |
543 | } | |
544 | } | |
545 | } | |
546 | } | |
7340864d | 547 | } |
70c95f95 | 548 | |
cd1c2af1 | 549 | //________________________________________________________ |
550 | void AliITSVertexer3D::MarkUsedClusters(){ | |
551 | // Mark clusters of tracklets used in vertex claulation | |
552 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ | |
553 | AliStrLine *lin = (AliStrLine*)fLines.At(i); | |
554 | Int_t idClu1=lin->GetIdPoint(0); | |
555 | Int_t idClu2=lin->GetIdPoint(1); | |
556 | fUsedCluster.SetBitNumber(idClu1); | |
557 | fUsedCluster.SetBitNumber(idClu2); | |
558 | } | |
559 | } | |
560 | //________________________________________________________ | |
561 | Int_t AliITSVertexer3D::RemoveTracklets(){ | |
562 | // Remove trackelts close to first found vertex | |
563 | Double_t vert[3]={fVert3D.GetXv(),fVert3D.GetYv(),fVert3D.GetZv()}; | |
564 | Int_t nRemoved=0; | |
565 | for(Int_t i=0; i<fLines.GetEntriesFast();i++){ | |
566 | AliStrLine *lin = (AliStrLine*)fLines.At(i); | |
567 | if(lin->GetDistFromPoint(vert)<fDCAforPileup){ | |
568 | Int_t idClu1=lin->GetIdPoint(0); | |
569 | Int_t idClu2=lin->GetIdPoint(1); | |
570 | fUsedCluster.SetBitNumber(idClu1); | |
571 | fUsedCluster.SetBitNumber(idClu2); | |
572 | fLines.RemoveAt(i); | |
573 | ++nRemoved; | |
574 | } | |
575 | } | |
576 | fLines.Compress(); | |
577 | return nRemoved; | |
578 | } | |
579 | //________________________________________________________ | |
580 | void AliITSVertexer3D::FindOther3DVertices(TTree *itsClusterTree){ | |
581 | // pileup identification based on 3D vertexing with not used clusters | |
582 | MarkUsedClusters(); | |
583 | fLines.Clear("C"); | |
584 | Int_t nolines = FindTracklets(itsClusterTree,2); | |
585 | if(nolines>=2){ | |
586 | Int_t nr=RemoveTracklets(); | |
587 | nolines-=nr; | |
588 | if(nolines>=2){ | |
589 | Int_t rc=Prepare3DVertex(2); | |
590 | if(rc==0){ | |
591 | fVert3D=AliVertexerTracks::TrackletVertexFinder(&fLines,0); | |
592 | if(fVert3D.GetNContributors()>=fMinTrackletsForPilup){ | |
593 | fIsPileup=kTRUE; | |
594 | fZpuv=fVert3D.GetZv(); | |
595 | fNTrpuv=fVert3D.GetNContributors(); | |
596 | } | |
597 | } | |
598 | } | |
599 | } | |
600 | } | |
601 | //______________________________________________________________________ | |
602 | void AliITSVertexer3D::PileupFromZ(){ | |
603 | // Calls the pileup algorithm of ALiITSVertexerZ | |
604 | Int_t binmin, binmax; | |
605 | Int_t nPeaks=AliITSVertexerZ::GetPeakRegion(fZHisto,binmin,binmax); | |
606 | if(nPeaks==2)AliWarning("2 peaks found"); | |
607 | Int_t firstPeakCont=0; | |
608 | Float_t firstPeakPos=0.; | |
609 | for(Int_t i=binmin-1;i<=binmax+1;i++){ | |
610 | firstPeakCont+=fZHisto->GetBinContent(i); | |
611 | firstPeakPos+=fZHisto->GetBinContent(i)*fZHisto->GetBinCenter(i); | |
612 | } | |
613 | firstPeakPos/=firstPeakCont; | |
614 | Int_t ncontr2=0; | |
615 | if(firstPeakCont>fMinTrackletsForPilup){ | |
616 | Float_t secPeakPos; | |
617 | ncontr2=AliITSVertexerZ::FindSecondPeak(fZHisto,binmin,binmax,secPeakPos); | |
618 | if(ncontr2>=fMinTrackletsForPilup){ | |
619 | fIsPileup=kTRUE; | |
620 | fZpuv=secPeakPos; | |
621 | fNTrpuv=ncontr2; | |
622 | } | |
623 | } | |
624 | } | |
70c95f95 | 625 | //________________________________________________________ |
626 | void AliITSVertexer3D::PrintStatus() const { | |
627 | // Print current status | |
6a0d56b8 | 628 | printf("=======================================================\n"); |
629 | printf("Loose cut on Delta Phi %f\n",fCoarseDiffPhiCut); | |
630 | printf("Cut on tracklet DCA to Z axis %f\n",fCoarseMaxRCut); | |
631 | printf("Cut on tracklet DCA to beam axis %f\n",fMaxRCut); | |
632 | printf("Cut on diamond (Z) %f\n",fZCutDiamond); | |
633 | printf("Cut on DCA - tracklet to tracklet and to vertex %f\n",fDCAcut); | |
cd1c2af1 | 634 | printf("Max Phi difference: %f\n",fDiffPhiMax); |
635 | printf("Pileup algo: %d\n",fPileupAlgo); | |
636 | printf("Min DCA to 1st vetrtex for pileup: %f\n",fDCAforPileup); | |
6a0d56b8 | 637 | printf("=======================================================\n"); |
70c95f95 | 638 | } |