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