]>
Commit | Line | Data |
---|---|---|
2d57349e | 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 | ||
16 | //----------------------------------------------------------------- | |
6d8df534 | 17 | // Implementation of the vertexer from tracks |
2d57349e | 18 | // |
19 | // Origin: AliITSVertexerTracks | |
20 | // A.Dainese, Padova, | |
21 | // andrea.dainese@pd.infn.it | |
22 | // M.Masera, Torino, | |
23 | // massimo.masera@to.infn.it | |
24 | // Moved to STEER and adapted to ESD tracks: | |
6d8df534 | 25 | // F.Prino, Torino, prino@to.infn.it |
2d57349e | 26 | //----------------------------------------------------------------- |
27 | ||
28 | //---- Root headers -------- | |
817e1004 | 29 | #include <TSystem.h> |
30 | #include <TClonesArray.h> | |
a7ede274 | 31 | #include <TDirectory.h> |
32 | #include <TFile.h> | |
2d57349e | 33 | //---- AliRoot headers ----- |
34 | #include "AliStrLine.h" | |
6d8df534 | 35 | #include "AliExternalTrackParam.h" |
2258e165 | 36 | #include "AliNeutralTrackParam.h" |
0aa2faf4 | 37 | #include "AliVEvent.h" |
38 | #include "AliVTrack.h" | |
2d57349e | 39 | #include "AliESDtrack.h" |
6d8df534 | 40 | #include "AliVertexerTracks.h" |
2d57349e | 41 | |
42 | ClassImp(AliVertexerTracks) | |
43 | ||
44 | ||
45 | //---------------------------------------------------------------------------- | |
46 | AliVertexerTracks::AliVertexerTracks(): | |
07680cae | 47 | TObject(), |
48 | fVert(), | |
49 | fCurrentVertex(0), | |
f09c879d | 50 | fMode(0), |
dc3719f3 | 51 | fFieldkG(-999.), |
6d8df534 | 52 | fTrkArraySel(), |
53 | fIdSel(0), | |
07680cae | 54 | fTrksToSkip(0), |
55 | fNTrksToSkip(0), | |
f09c879d | 56 | fConstraint(kFALSE), |
57 | fOnlyFitter(kFALSE), | |
58 | fMinTracks(1), | |
6b29d399 | 59 | fMinClusters(5), |
f09c879d | 60 | fDCAcut(0.1), |
61 | fDCAcutIter0(0.1), | |
62 | fNSigma(3.), | |
fee9ef0d | 63 | fMaxd0z0(0.5), |
f09c879d | 64 | fMinDetFitter(100.), |
65 | fMaxTgl(1000.), | |
50ff0bcd | 66 | fITSrefit(kTRUE), |
40cfa14d | 67 | fITSpureSA(kFALSE), |
6b29d399 | 68 | fFiducialR(3.), |
69 | fFiducialZ(30.), | |
3b768a1d | 70 | fnSigmaForUi00(1.5), |
8c75f668 | 71 | fAlgo(1), |
ce1c94a7 | 72 | fAlgoIter0(4), |
73 | fSelectOnTOFBunchCrossing(kFALSE), | |
74 | fKeepAlsoUnflaggedTOFBunchCrossing(kTRUE) | |
2d57349e | 75 | { |
76 | // | |
77 | // Default constructor | |
78 | // | |
3b768a1d | 79 | SetVtxStart(); |
07680cae | 80 | SetVtxStartSigma(); |
2d57349e | 81 | } |
dc3719f3 | 82 | //---------------------------------------------------------------------------- |
83 | AliVertexerTracks::AliVertexerTracks(Double_t fieldkG): | |
07680cae | 84 | TObject(), |
85 | fVert(), | |
86 | fCurrentVertex(0), | |
f09c879d | 87 | fMode(0), |
88 | fFieldkG(fieldkG), | |
6d8df534 | 89 | fTrkArraySel(), |
90 | fIdSel(0), | |
07680cae | 91 | fTrksToSkip(0), |
92 | fNTrksToSkip(0), | |
f09c879d | 93 | fConstraint(kFALSE), |
94 | fOnlyFitter(kFALSE), | |
95 | fMinTracks(1), | |
6b29d399 | 96 | fMinClusters(5), |
f09c879d | 97 | fDCAcut(0.1), |
98 | fDCAcutIter0(0.1), | |
99 | fNSigma(3.), | |
fee9ef0d | 100 | fMaxd0z0(0.5), |
f09c879d | 101 | fMinDetFitter(100.), |
102 | fMaxTgl(1000.), | |
50ff0bcd | 103 | fITSrefit(kTRUE), |
40cfa14d | 104 | fITSpureSA(kFALSE), |
6b29d399 | 105 | fFiducialR(3.), |
106 | fFiducialZ(30.), | |
3b768a1d | 107 | fnSigmaForUi00(1.5), |
8c75f668 | 108 | fAlgo(1), |
ce1c94a7 | 109 | fAlgoIter0(4), |
110 | fSelectOnTOFBunchCrossing(kFALSE), | |
111 | fKeepAlsoUnflaggedTOFBunchCrossing(kTRUE) | |
2d57349e | 112 | { |
113 | // | |
114 | // Standard constructor | |
115 | // | |
3b768a1d | 116 | SetVtxStart(); |
07680cae | 117 | SetVtxStartSigma(); |
2d57349e | 118 | } |
119 | //----------------------------------------------------------------------------- | |
fee9ef0d | 120 | AliVertexerTracks::~AliVertexerTracks() |
121 | { | |
2d57349e | 122 | // Default Destructor |
eae1677e | 123 | // The objects pointed by the following pointer are not owned |
2d57349e | 124 | // by this class and are not deleted |
ec1be5d5 | 125 | fCurrentVertex = 0; |
6fefa5ce | 126 | if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; } |
127 | if(fIdSel) { delete [] fIdSel; fIdSel=NULL; } | |
ec1be5d5 | 128 | } |
c5e3e5d1 | 129 | //---------------------------------------------------------------------------- |
f5740e9a | 130 | AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent) |
ec1be5d5 | 131 | { |
c5e3e5d1 | 132 | // |
0aa2faf4 | 133 | // Primary vertex for current ESD or AOD event |
eae1677e | 134 | // (Two iterations: |
fee9ef0d | 135 | // 1st with 5*fNSigma*sigma cut w.r.t. to initial vertex |
136 | // + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE | |
137 | // 2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration) | |
07680cae | 138 | // |
139 | fCurrentVertex = 0; | |
140 | ||
0aa2faf4 | 141 | TString evtype = vEvent->IsA()->GetName(); |
142 | Bool_t inputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE); | |
143 | ||
144 | if(inputAOD && fMode==1) { | |
145 | printf("Error : AliVertexerTracks: no TPC-only vertex from AOD\n"); | |
146 | TooFewTracks(); | |
147 | return fCurrentVertex; | |
148 | } | |
149 | ||
dc3719f3 | 150 | // accept 1-track case only if constraint is available |
151 | if(!fConstraint && fMinTracks==1) fMinTracks=2; | |
152 | ||
0aa2faf4 | 153 | // read tracks from AlivEvent |
154 | Int_t nTrks = (Int_t)vEvent->GetNumberOfTracks(); | |
f09c879d | 155 | if(nTrks<fMinTracks) { |
6d8df534 | 156 | TooFewTracks(); |
80d36431 | 157 | return fCurrentVertex; |
158 | } | |
6d8df534 | 159 | // |
713e5fbf | 160 | |
a7ede274 | 161 | TDirectory * olddir = gDirectory; |
f09c879d | 162 | TFile *f = 0; |
163 | if(nTrks>500) f = new TFile("VertexerTracks.root","recreate"); | |
6d8df534 | 164 | TObjArray trkArrayOrig(nTrks); |
165 | UShort_t *idOrig = new UShort_t[nTrks]; | |
07680cae | 166 | |
6d8df534 | 167 | Int_t nTrksOrig=0; |
f09c879d | 168 | AliExternalTrackParam *t=0; |
0aa2faf4 | 169 | // loop on tracks |
6d8df534 | 170 | for(Int_t i=0; i<nTrks; i++) { |
0aa2faf4 | 171 | AliVTrack *track = (AliVTrack*)vEvent->GetTrack(i); |
4539449a | 172 | // check tracks to skip |
173 | Bool_t skipThis = kFALSE; | |
174 | for(Int_t j=0; j<fNTrksToSkip; j++) { | |
175 | if(track->GetID()==fTrksToSkip[j]) { | |
176 | AliDebug(1,Form("skipping track: %d",i)); | |
177 | skipThis = kTRUE; | |
178 | } | |
179 | } | |
180 | if(skipThis) continue; | |
0aa2faf4 | 181 | |
40cfa14d | 182 | // skip pure ITS SA tracks (default) |
183 | if(!fITSpureSA && (track->GetStatus()&AliESDtrack::kITSpureSA)) continue; | |
184 | // or use only pure ITS SA tracks | |
185 | if(fITSpureSA && !(track->GetStatus()&AliESDtrack::kITSpureSA)) continue; | |
71d7871c | 186 | |
02bb7b6c | 187 | // kITSrefit |
188 | if(fMode==0 && fITSrefit && !(track->GetStatus()&AliESDtrack::kITSrefit)) continue; | |
189 | ||
d8e3340c | 190 | // reject tracks according to bunch crossing id from TOF |
191 | if(fSelectOnTOFBunchCrossing) { | |
192 | Int_t bcTOF = track->GetTOFBunchCrossing(); | |
193 | if(bcTOF>0) continue; | |
194 | if(bcTOF==-1 && !fKeepAlsoUnflaggedTOFBunchCrossing) continue; | |
195 | } | |
196 | ||
197 | ||
4539449a | 198 | if(!inputAOD) { // ESD |
0aa2faf4 | 199 | AliESDtrack* esdt = (AliESDtrack*)track; |
200 | if(esdt->GetNcls(fMode) < fMinClusters) continue; | |
4539449a | 201 | if(fMode==0) { // ITS mode |
4539449a | 202 | Double_t x,p[5],cov[15]; |
203 | esdt->GetExternalParameters(x,p); | |
204 | esdt->GetExternalCovariance(cov); | |
205 | t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov); | |
206 | } else if(fMode==1) { // TPC mode | |
207 | t = (AliExternalTrackParam*)esdt->GetTPCInnerParam(); | |
208 | if(!t) continue; | |
209 | Double_t radius = 2.8; //something less than the beam pipe radius | |
210 | if(!PropagateTrackTo(t,radius)) continue; | |
211 | } | |
212 | } else { // AOD (only ITS mode) | |
213 | Int_t ncls0=0; | |
214 | for(Int_t l=0;l<6;l++) if(TESTBIT(track->GetITSClusterMap(),l)) ncls0++; | |
215 | if(ncls0 < fMinClusters) continue; | |
0aa2faf4 | 216 | t = new AliExternalTrackParam(track); |
50ff0bcd | 217 | } |
6d8df534 | 218 | trkArrayOrig.AddLast(t); |
0aa2faf4 | 219 | idOrig[nTrksOrig]=(UShort_t)track->GetID(); |
6d8df534 | 220 | nTrksOrig++; |
0aa2faf4 | 221 | } // end loop on tracks |
fee9ef0d | 222 | |
f09c879d | 223 | // call method that will reconstruct the vertex |
6d8df534 | 224 | FindPrimaryVertex(&trkArrayOrig,idOrig); |
225 | ||
4539449a | 226 | if(fMode==0) trkArrayOrig.Delete(); |
6fefa5ce | 227 | delete [] idOrig; idOrig=NULL; |
f09c879d | 228 | |
229 | if(f) { | |
230 | f->Close(); delete f; f = NULL; | |
231 | gSystem->Unlink("VertexerTracks.root"); | |
232 | olddir->cd(); | |
233 | } | |
6d8df534 | 234 | |
4448efe4 | 235 | // set vertex ID for tracks used in the fit |
236 | // (only for ESD) | |
237 | if(!inputAOD) { | |
238 | Int_t nIndices = fCurrentVertex->GetNIndices(); | |
239 | UShort_t *indices = fCurrentVertex->GetIndices(); | |
240 | for(Int_t ind=0; ind<nIndices; ind++) { | |
241 | AliESDtrack *esdt = (AliESDtrack*)vEvent->GetTrack(indices[ind]); | |
242 | esdt->SetVertexID(-1); | |
243 | } | |
244 | } | |
245 | ||
6d8df534 | 246 | return fCurrentVertex; |
247 | } | |
248 | //---------------------------------------------------------------------------- | |
f5740e9a | 249 | AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const TObjArray *trkArrayOrig, |
6d8df534 | 250 | UShort_t *idOrig) |
251 | { | |
252 | // | |
253 | // Primary vertex using the AliExternalTrackParam's in the TObjArray. | |
254 | // idOrig must contain the track IDs from AliESDtrack::GetID() | |
255 | // (Two iterations: | |
256 | // 1st with 5*fNSigma*sigma cut w.r.t. to initial vertex | |
257 | // + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE | |
258 | // 2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration) | |
259 | // | |
260 | fCurrentVertex = 0; | |
261 | ||
262 | // accept 1-track case only if constraint is available | |
263 | if(!fConstraint && fMinTracks==1) fMinTracks=2; | |
264 | ||
f09c879d | 265 | // read tracks from array |
6d8df534 | 266 | Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast(); |
919e537f | 267 | AliDebug(1,Form("Initial number of tracks: %d",nTrksOrig)); |
f09c879d | 268 | if(nTrksOrig<fMinTracks) { |
919e537f | 269 | AliDebug(1,"TooFewTracks"); |
6d8df534 | 270 | TooFewTracks(); |
271 | return fCurrentVertex; | |
272 | } | |
273 | ||
fee9ef0d | 274 | // If fConstraint=kFALSE |
275 | // run VertexFinder(1) to get rough estimate of initVertex (x,y) | |
276 | if(!fConstraint) { | |
6d8df534 | 277 | // fill fTrkArraySel, for VertexFinder() |
278 | fIdSel = new UShort_t[nTrksOrig]; | |
279 | PrepareTracks(*trkArrayOrig,idOrig,0); | |
6fefa5ce | 280 | if(fIdSel) { delete [] fIdSel; fIdSel=NULL; } |
6d8df534 | 281 | Double_t cutsave = fDCAcut; |
f09c879d | 282 | fDCAcut = fDCAcutIter0; |
8c75f668 | 283 | // vertex finder |
284 | switch (fAlgoIter0) { | |
285 | case 1: StrLinVertexFinderMinDist(1); break; | |
286 | case 2: StrLinVertexFinderMinDist(0); break; | |
287 | case 3: HelixVertexFinder(); break; | |
288 | case 4: VertexFinder(1); break; | |
289 | case 5: VertexFinder(0); break; | |
290 | default: printf("Wrong algorithm\n"); break; | |
291 | } | |
fee9ef0d | 292 | fDCAcut = cutsave; |
fee9ef0d | 293 | if(fVert.GetNContributors()>0) { |
294 | fVert.GetXYZ(fNominalPos); | |
295 | fNominalPos[0] = fVert.GetXv(); | |
296 | fNominalPos[1] = fVert.GetYv(); | |
297 | fNominalPos[2] = fVert.GetZv(); | |
919e537f | 298 | AliDebug(1,Form("No mean vertex: VertexFinder gives (%f, %f, %f)",fNominalPos[0],fNominalPos[1],fNominalPos[2])); |
fee9ef0d | 299 | } else { |
300 | fNominalPos[0] = 0.; | |
301 | fNominalPos[1] = 0.; | |
302 | fNominalPos[2] = 0.; | |
919e537f | 303 | AliDebug(1,"No mean vertex and VertexFinder failed"); |
0c69be49 | 304 | } |
07680cae | 305 | } |
fee9ef0d | 306 | |
307 | // TWO ITERATIONS: | |
308 | // | |
309 | // ITERATION 1 | |
310 | // propagate tracks to fNominalPos vertex | |
311 | // preselect them: | |
312 | // if(constraint) reject for |d0|>5*fNSigma*sigma w.r.t. fNominal... vertex | |
313 | // else reject for |d0|\oplus|z0| > 5 mm w.r.t. fNominal... vertex | |
07680cae | 314 | // ITERATION 2 |
315 | // propagate tracks to best between initVertex and fCurrentVertex | |
316 | // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best | |
317 | // between initVertex and fCurrentVertex) | |
6d8df534 | 318 | for(Int_t iter=1; iter<=2; iter++) { |
319 | if(fOnlyFitter && iter==1) continue; | |
6fefa5ce | 320 | if(fIdSel) { delete [] fIdSel; fIdSel=NULL; } |
6d8df534 | 321 | fIdSel = new UShort_t[nTrksOrig]; |
322 | Int_t nTrksSel = PrepareTracks(*trkArrayOrig,idOrig,iter); | |
919e537f | 323 | AliDebug(1,Form("N tracks selected in iteration %d: %d",iter,nTrksSel)); |
6d8df534 | 324 | if(nTrksSel < fMinTracks) { |
325 | TooFewTracks(); | |
fee9ef0d | 326 | return fCurrentVertex; |
327 | } | |
07680cae | 328 | |
fee9ef0d | 329 | // vertex finder |
330 | if(!fOnlyFitter) { | |
6d8df534 | 331 | if(nTrksSel==1) { |
919e537f | 332 | AliDebug(1,"Just one track"); |
fee9ef0d | 333 | OneTrackVertFinder(); |
6d8df534 | 334 | } else { |
fee9ef0d | 335 | switch (fAlgo) { |
336 | case 1: StrLinVertexFinderMinDist(1); break; | |
337 | case 2: StrLinVertexFinderMinDist(0); break; | |
338 | case 3: HelixVertexFinder(); break; | |
339 | case 4: VertexFinder(1); break; | |
340 | case 5: VertexFinder(0); break; | |
919e537f | 341 | default: printf("Wrong algorithm"); break; |
fee9ef0d | 342 | } |
343 | } | |
919e537f | 344 | AliDebug(1," Vertex finding completed"); |
0c69be49 | 345 | } |
07680cae | 346 | |
fee9ef0d | 347 | // vertex fitter |
6d8df534 | 348 | VertexFitter(); |
fee9ef0d | 349 | } // end loop on the two iterations |
07680cae | 350 | |
fee9ef0d | 351 | |
50ff0bcd | 352 | // set indices of used tracks |
353 | UShort_t *indices = 0; | |
50ff0bcd | 354 | if(fCurrentVertex->GetNContributors()>0) { |
4d023a8f | 355 | Int_t nIndices = (Int_t)fTrkArraySel.GetEntriesFast(); |
356 | indices = new UShort_t[nIndices]; | |
357 | for(Int_t jj=0; jj<nIndices; jj++) | |
6d8df534 | 358 | indices[jj] = fIdSel[jj]; |
4d023a8f | 359 | fCurrentVertex->SetIndices(nIndices,indices); |
50ff0bcd | 360 | } |
6fefa5ce | 361 | if (indices) {delete [] indices; indices=NULL;} |
6d8df534 | 362 | // |
fee9ef0d | 363 | |
6d8df534 | 364 | // set vertex title |
365 | TString title="VertexerTracksNoConstraint"; | |
366 | if(fConstraint) { | |
367 | title="VertexerTracksWithConstraint"; | |
368 | if(fOnlyFitter) title.Append("OnlyFitter"); | |
369 | } | |
370 | fCurrentVertex->SetTitle(title.Data()); | |
371 | // | |
a7ede274 | 372 | |
919e537f | 373 | AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors())); |
50ff0bcd | 374 | |
6d8df534 | 375 | // clean up |
6fefa5ce | 376 | delete [] fIdSel; fIdSel=NULL; |
6d8df534 | 377 | fTrkArraySel.Delete(); |
6fefa5ce | 378 | if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; } |
6d8df534 | 379 | // |
fee9ef0d | 380 | |
07680cae | 381 | return fCurrentVertex; |
382 | } | |
50ff0bcd | 383 | //------------------------------------------------------------------------ |
fee9ef0d | 384 | Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3]) |
385 | { | |
50ff0bcd | 386 | // |
387 | Double_t det=matr[0][0]*matr[1][1]*matr[2][2]-matr[0][0]*matr[1][2]*matr[2][1]-matr[0][1]*matr[1][0]*matr[2][2]+matr[0][1]*matr[1][2]*matr[2][0]+matr[0][2]*matr[1][0]*matr[2][1]-matr[0][2]*matr[1][1]*matr[2][0]; | |
388 | return det; | |
389 | } | |
390 | //------------------------------------------------------------------------- | |
f5740e9a | 391 | void AliVertexerTracks::GetStrLinDerivMatrix(const Double_t *p0,const Double_t *p1,Double_t (*m)[3],Double_t *d) |
fee9ef0d | 392 | { |
50ff0bcd | 393 | // |
394 | Double_t x12=p0[0]-p1[0]; | |
395 | Double_t y12=p0[1]-p1[1]; | |
396 | Double_t z12=p0[2]-p1[2]; | |
397 | Double_t kk=x12*x12+y12*y12+z12*z12; | |
398 | m[0][0]=2-2/kk*x12*x12; | |
399 | m[0][1]=-2/kk*x12*y12; | |
400 | m[0][2]=-2/kk*x12*z12; | |
401 | m[1][0]=-2/kk*x12*y12; | |
402 | m[1][1]=2-2/kk*y12*y12; | |
403 | m[1][2]=-2/kk*y12*z12; | |
404 | m[2][0]=-2/kk*x12*z12; | |
405 | m[2][1]=-2*y12*z12; | |
406 | m[2][2]=2-2/kk*z12*z12; | |
407 | d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12; | |
408 | d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12; | |
409 | d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12; | |
c5e3e5d1 | 410 | |
50ff0bcd | 411 | } |
412 | //-------------------------------------------------------------------------- | |
f5740e9a | 413 | void AliVertexerTracks::GetStrLinDerivMatrix(const Double_t *p0,const Double_t *p1,const Double_t *sigmasq,Double_t (*m)[3],Double_t *d) |
fee9ef0d | 414 | { |
50ff0bcd | 415 | // |
416 | Double_t x12=p1[0]-p0[0]; | |
417 | Double_t y12=p1[1]-p0[1]; | |
418 | Double_t z12=p1[2]-p0[2]; | |
c5e3e5d1 | 419 | |
50ff0bcd | 420 | Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1]; |
c5e3e5d1 | 421 | |
50ff0bcd | 422 | Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]); |
c5e3e5d1 | 423 | |
50ff0bcd | 424 | Double_t cc[3]; |
425 | cc[0]=-x12/sigmasq[0]; | |
426 | cc[1]=-y12/sigmasq[1]; | |
427 | cc[2]=-z12/sigmasq[2]; | |
ec1be5d5 | 428 | |
50ff0bcd | 429 | Double_t ww=(-p0[0]*x12*sigmasq[1]*sigmasq[2]-p0[1]*y12*sigmasq[0]*sigmasq[2]-p0[2]*z12*sigmasq[0]*sigmasq[1])/den; |
ec1be5d5 | 430 | |
50ff0bcd | 431 | Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2]; |
2d57349e | 432 | |
50ff0bcd | 433 | Double_t aa[3]; |
434 | aa[0]=x12*sigmasq[1]*sigmasq[2]/den; | |
435 | aa[1]=y12*sigmasq[0]*sigmasq[2]/den; | |
436 | aa[2]=z12*sigmasq[0]*sigmasq[1]/den; | |
2d57349e | 437 | |
50ff0bcd | 438 | m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0]; |
439 | m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0]; | |
440 | m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0]; | |
2d57349e | 441 | |
50ff0bcd | 442 | m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1]; |
443 | m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1]; | |
444 | m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1]; | |
2d57349e | 445 | |
50ff0bcd | 446 | m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2]; |
447 | m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2]; | |
448 | m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2]; | |
08df6187 | 449 | |
50ff0bcd | 450 | d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0]; |
451 | d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1]; | |
452 | d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2]; | |
2d57349e | 453 | |
2d57349e | 454 | } |
50ff0bcd | 455 | //-------------------------------------------------------------------------- |
f5740e9a | 456 | Double_t AliVertexerTracks::GetStrLinMinDist(const Double_t *p0,const Double_t *p1,const Double_t *x0) |
fee9ef0d | 457 | { |
50ff0bcd | 458 | // |
459 | Double_t x12=p0[0]-p1[0]; | |
460 | Double_t y12=p0[1]-p1[1]; | |
461 | Double_t z12=p0[2]-p1[2]; | |
462 | Double_t x10=p0[0]-x0[0]; | |
463 | Double_t y10=p0[1]-x0[1]; | |
464 | Double_t z10=p0[2]-x0[2]; | |
3446476f | 465 | // return ((x10*x10+y10*y10+z10*z10)*(x12*x12+y12*y12+z12*z12)-(x10*x12+y10*y12+z10*z12)*(x10*x12+y10*y12+z10*z12))/(x12*x12+y12*y12+z12*z12); |
3446476f | 466 | return ((y10*z12-z10*y12)*(y10*z12-z10*y12)+ |
467 | (z10*x12-x10*z12)*(z10*x12-x10*z12)+ | |
468 | (x10*y12-y10*x12)*(x10*y12-y10*x12)) | |
469 | /(x12*x12+y12*y12+z12*z12); | |
2d57349e | 470 | } |
471 | //--------------------------------------------------------------------------- | |
fee9ef0d | 472 | void AliVertexerTracks::OneTrackVertFinder() |
473 | { | |
0c69be49 | 474 | // find vertex for events with 1 track, using DCA to nominal beam axis |
919e537f | 475 | AliDebug(1,Form("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries())); |
6d8df534 | 476 | AliExternalTrackParam *track1; |
477 | track1 = (AliExternalTrackParam*)fTrkArraySel.At(0); | |
0c69be49 | 478 | Double_t alpha=track1->GetAlpha(); |
479 | Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1]; | |
480 | Double_t pos[3],dir[3]; | |
f09c879d | 481 | track1->GetXYZAt(mindist,GetFieldkG(),pos); |
482 | track1->GetPxPyPzAt(mindist,GetFieldkG(),dir); | |
0c69be49 | 483 | AliStrLine *line1 = new AliStrLine(pos,dir); |
484 | Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.}; | |
485 | Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.}; | |
486 | AliStrLine *zeta=new AliStrLine(p1,p2,kTRUE); | |
487 | Double_t crosspoint[3]={0.,0.,0.}; | |
488 | Double_t sigma=999.; | |
489 | Int_t nContrib=-1; | |
490 | Int_t retcode = zeta->Cross(line1,crosspoint); | |
491 | if(retcode>=0){ | |
492 | sigma=line1->GetDistFromPoint(crosspoint); | |
493 | nContrib=1; | |
494 | } | |
495 | delete zeta; | |
496 | delete line1; | |
497 | fVert.SetXYZ(crosspoint); | |
498 | fVert.SetDispersion(sigma); | |
499 | fVert.SetNContributors(nContrib); | |
500 | } | |
501 | //--------------------------------------------------------------------------- | |
fee9ef0d | 502 | void AliVertexerTracks::HelixVertexFinder() |
503 | { | |
2d57349e | 504 | // Get estimate of vertex position in (x,y) from tracks DCA |
505 | ||
506 | ||
507 | Double_t initPos[3]; | |
508 | initPos[2] = 0.; | |
509 | for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i]; | |
2d57349e | 510 | |
6d8df534 | 511 | Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast(); |
2d57349e | 512 | |
513 | Double_t aver[3]={0.,0.,0.}; | |
514 | Double_t averquad[3]={0.,0.,0.}; | |
515 | Double_t sigmaquad[3]={0.,0.,0.}; | |
516 | Double_t sigma=0; | |
517 | Int_t ncombi = 0; | |
6d8df534 | 518 | AliExternalTrackParam *track1; |
519 | AliExternalTrackParam *track2; | |
2d57349e | 520 | Double_t distCA; |
6d8df534 | 521 | Double_t x; |
2d57349e | 522 | Double_t alpha, cs, sn; |
523 | Double_t crosspoint[3]; | |
524 | for(Int_t i=0; i<nacc; i++){ | |
6d8df534 | 525 | track1 = (AliExternalTrackParam*)fTrkArraySel.At(i); |
2d57349e | 526 | |
527 | ||
528 | for(Int_t j=i+1; j<nacc; j++){ | |
6d8df534 | 529 | track2 = (AliExternalTrackParam*)fTrkArraySel.At(j); |
2d57349e | 530 | |
f09c879d | 531 | distCA=track2->PropagateToDCA(track1,GetFieldkG()); |
2d57349e | 532 | if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){ |
6d8df534 | 533 | x=track1->GetX(); |
2d57349e | 534 | alpha=track1->GetAlpha(); |
535 | cs=TMath::Cos(alpha); sn=TMath::Sin(alpha); | |
6d8df534 | 536 | Double_t x1=x*cs - track1->GetY()*sn; |
537 | Double_t y1=x*sn + track1->GetY()*cs; | |
538 | Double_t z1=track1->GetZ(); | |
539 | ||
2d57349e | 540 | Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2(); |
6d8df534 | 541 | x=track2->GetX(); |
2d57349e | 542 | alpha=track2->GetAlpha(); |
543 | cs=TMath::Cos(alpha); sn=TMath::Sin(alpha); | |
6d8df534 | 544 | Double_t x2=x*cs - track2->GetY()*sn; |
545 | Double_t y2=x*sn + track2->GetY()*cs; | |
546 | Double_t z2=track2->GetZ(); | |
2d57349e | 547 | Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2(); |
548 | Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2(); | |
549 | Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1; | |
550 | Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1; | |
551 | Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1; | |
552 | crosspoint[0]=wx1*x1 + wx2*x2; | |
553 | crosspoint[1]=wy1*y1 + wy2*y2; | |
554 | crosspoint[2]=wz1*z1 + wz2*z2; | |
555 | ||
556 | ncombi++; | |
557 | for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj]; | |
558 | for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]); | |
559 | } | |
560 | } | |
561 | ||
562 | } | |
563 | if(ncombi>0){ | |
564 | for(Int_t jj=0;jj<3;jj++){ | |
565 | initPos[jj] = aver[jj]/ncombi; | |
566 | averquad[jj]/=ncombi; | |
567 | sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj]; | |
568 | sigma+=sigmaquad[jj]; | |
569 | } | |
570 | sigma=TMath::Sqrt(TMath::Abs(sigma)); | |
571 | } | |
572 | else { | |
573 | Warning("HelixVertexFinder","Finder did not succed"); | |
574 | sigma=999; | |
575 | } | |
576 | fVert.SetXYZ(initPos); | |
577 | fVert.SetDispersion(sigma); | |
578 | fVert.SetNContributors(ncombi); | |
579 | } | |
50ff0bcd | 580 | //---------------------------------------------------------------------------- |
f5740e9a | 581 | Int_t AliVertexerTracks::PrepareTracks(const TObjArray &trkArrayOrig, |
582 | const UShort_t *idOrig, | |
6d8df534 | 583 | Int_t optImpParCut) |
fee9ef0d | 584 | { |
50ff0bcd | 585 | // |
586 | // Propagate tracks to initial vertex position and store them in a TObjArray | |
587 | // | |
919e537f | 588 | AliDebug(1," PrepareTracks()"); |
6d8df534 | 589 | |
590 | Int_t nTrksOrig = (Int_t)trkArrayOrig.GetEntriesFast(); | |
591 | Int_t nTrksSel = 0; | |
fee9ef0d | 592 | Double_t maxd0rphi; |
50ff0bcd | 593 | Double_t sigmaCurr[3]; |
fee9ef0d | 594 | Double_t normdistx,normdisty; |
6d8df534 | 595 | Double_t d0z0[2],covd0z0[3]; |
43c9dae1 | 596 | Double_t sigmad0; |
50ff0bcd | 597 | |
598 | AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1); | |
599 | ||
6d8df534 | 600 | if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete(); |
50ff0bcd | 601 | |
6d8df534 | 602 | AliExternalTrackParam *track=0; |
50ff0bcd | 603 | |
f09c879d | 604 | // loop on tracks |
6d8df534 | 605 | for(Int_t i=0; i<nTrksOrig; i++) { |
2258e165 | 606 | AliExternalTrackParam *trackOrig=(AliExternalTrackParam*)trkArrayOrig.At(i); |
607 | if(trackOrig->Charge()!=0) { // normal tracks | |
608 | track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i)); | |
609 | } else { // neutral tracks (from a V0) | |
610 | track = new AliNeutralTrackParam(*(AliNeutralTrackParam*)trkArrayOrig.At(i)); | |
611 | } | |
02bb7b6c | 612 | |
f09c879d | 613 | // tgl cut |
614 | if(TMath::Abs(track->GetTgl())>fMaxTgl) { | |
919e537f | 615 | AliDebug(1,Form(" rejecting track with tgl = %f",track->GetTgl())); |
f09c879d | 616 | delete track; continue; |
6d8df534 | 617 | } |
50ff0bcd | 618 | |
d993e038 | 619 | Bool_t propagateOK = kFALSE, cutond0z0 = kTRUE; |
50ff0bcd | 620 | // propagate track to vertex |
43c9dae1 | 621 | if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0 |
18eba197 | 622 | propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0); |
fee9ef0d | 623 | } else { // optImpParCut==2 |
50ff0bcd | 624 | fCurrentVertex->GetSigmaXYZ(sigmaCurr); |
fee9ef0d | 625 | normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]); |
80d36431 | 626 | normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]); |
d993e038 | 627 | AliDebug(1,Form("normdistx %f %f %f",fCurrentVertex->GetXv(),fNominalPos[0],TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]))); |
628 | AliDebug(1,Form("normdisty %f %f %f",fCurrentVertex->GetYv(),fNominalPos[1],TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]))); | |
629 | AliDebug(1,Form("sigmaCurr %f %f %f",sigmaCurr[0],sigmaCurr[1],TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))); | |
fee9ef0d | 630 | if(normdistx < 3. && normdisty < 3. && |
80d36431 | 631 | (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) { |
18eba197 | 632 | propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0); |
50ff0bcd | 633 | } else { |
18eba197 | 634 | propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0); |
d993e038 | 635 | if(fConstraint) cutond0z0=kFALSE; |
50ff0bcd | 636 | } |
637 | } | |
638 | ||
18eba197 | 639 | if(!propagateOK) { |
919e537f | 640 | AliDebug(1," rejected"); |
18eba197 | 641 | delete track; continue; |
642 | } | |
643 | ||
43c9dae1 | 644 | sigmad0 = TMath::Sqrt(covd0z0[0]); |
645 | maxd0rphi = fNSigma*sigmad0; | |
fee9ef0d | 646 | if(optImpParCut==1) maxd0rphi *= 5.; |
6b29d399 | 647 | maxd0rphi = TMath::Min(maxd0rphi,fFiducialR); |
d993e038 | 648 | //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]); |
fee9ef0d | 649 | |
919e537f | 650 | AliDebug(1,Form("trk %d; id %d; |d0| = %f; d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f",i,(Int_t)idOrig[i],TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),fMaxd0z0)); |
f09c879d | 651 | |
652 | ||
653 | //---- track selection based on impact parameters ----// | |
654 | ||
655 | // always reject tracks outside fiducial volume | |
6b29d399 | 656 | if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) { |
919e537f | 657 | AliDebug(1," rejected"); |
f09c879d | 658 | delete track; continue; |
659 | } | |
660 | ||
661 | // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi | |
662 | if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) { | |
919e537f | 663 | AliDebug(1," rejected"); |
f09c879d | 664 | delete track; continue; |
665 | } | |
fee9ef0d | 666 | |
43c9dae1 | 667 | // if fConstraint=kFALSE, during iterations 1 and 2 || |
668 | // if fConstraint=kTRUE, during iteration 2, | |
f09c879d | 669 | // select tracks with d0oplusz0 < fMaxd0z0 |
43c9dae1 | 670 | if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) || |
d993e038 | 671 | ( fConstraint && optImpParCut==2 && cutond0z0)) { |
43c9dae1 | 672 | if(nTrksOrig>=3 && |
f09c879d | 673 | TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) { |
919e537f | 674 | AliDebug(1," rejected"); |
fee9ef0d | 675 | delete track; continue; |
676 | } | |
677 | } | |
43c9dae1 | 678 | |
f09c879d | 679 | // track passed all selections |
6d8df534 | 680 | fTrkArraySel.AddLast(track); |
681 | fIdSel[nTrksSel] = idOrig[i]; | |
682 | nTrksSel++; | |
f09c879d | 683 | } // end loop on tracks |
50ff0bcd | 684 | |
685 | delete initVertex; | |
686 | ||
6d8df534 | 687 | return nTrksSel; |
50ff0bcd | 688 | } |
f09c879d | 689 | //---------------------------------------------------------------------------- |
690 | Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track, | |
691 | Double_t xToGo) { | |
692 | //---------------------------------------------------------------- | |
693 | // COPIED from AliTracker | |
694 | // | |
695 | // Propagates the track to the plane X=xk (cm). | |
696 | // | |
697 | // Origin: Marian Ivanov, Marian.Ivanov@cern.ch | |
698 | //---------------------------------------------------------------- | |
699 | ||
700 | const Double_t kEpsilon = 0.00001; | |
701 | Double_t xpos = track->GetX(); | |
702 | Double_t dir = (xpos<xToGo) ? 1. : -1.; | |
703 | Double_t maxStep = 5; | |
704 | Double_t maxSnp = 0.8; | |
705 | // | |
706 | while ( (xToGo-xpos)*dir > kEpsilon){ | |
707 | Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep); | |
708 | Double_t x = xpos+step; | |
709 | Double_t xyz0[3],xyz1[3]; | |
710 | track->GetXYZ(xyz0); //starting global position | |
711 | ||
712 | if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE; // no prolongation | |
713 | xyz1[2]+=kEpsilon; // waiting for bug correction in geo | |
714 | ||
715 | if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE; | |
716 | if(!track->PropagateTo(x,GetFieldkG())) return kFALSE; | |
717 | ||
718 | if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE; | |
719 | track->GetXYZ(xyz0); // global position | |
720 | Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]); | |
721 | // | |
722 | Double_t ca=TMath::Cos(alphan-track->GetAlpha()), | |
723 | sa=TMath::Sin(alphan-track->GetAlpha()); | |
60e55aee | 724 | Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf)); |
f09c879d | 725 | Double_t sinNew = sf*ca - cf*sa; |
726 | if(TMath::Abs(sinNew) >= maxSnp) return kFALSE; | |
727 | if(!track->Rotate(alphan)) return kFALSE; | |
728 | ||
729 | xpos = track->GetX(); | |
730 | } | |
731 | return kTRUE; | |
732 | } | |
50ff0bcd | 733 | //--------------------------------------------------------------------------- |
fee9ef0d | 734 | AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx, |
f5740e9a | 735 | const TObjArray *trkArray, |
6d8df534 | 736 | UShort_t *id, |
f5740e9a | 737 | const Float_t *diamondxy) const |
fee9ef0d | 738 | { |
50ff0bcd | 739 | // |
fee9ef0d | 740 | // Removes tracks in trksTree from fit of inVtx |
50ff0bcd | 741 | // |
fee9ef0d | 742 | |
146c29df | 743 | if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) { |
919e537f | 744 | printf("ERROR: primary vertex has no constraint: cannot remove tracks"); |
146c29df | 745 | return 0x0; |
746 | } | |
fee9ef0d | 747 | if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter")) |
919e537f | 748 | printf("WARNING: result of tracks' removal will be only approximately correct"); |
fee9ef0d | 749 | |
750 | TMatrixD rv(3,1); | |
751 | rv(0,0) = inVtx->GetXv(); | |
752 | rv(1,0) = inVtx->GetYv(); | |
753 | rv(2,0) = inVtx->GetZv(); | |
754 | TMatrixD vV(3,3); | |
755 | Double_t cov[6]; | |
756 | inVtx->GetCovMatrix(cov); | |
757 | vV(0,0) = cov[0]; | |
758 | vV(0,1) = cov[1]; vV(1,0) = cov[1]; | |
759 | vV(1,1) = cov[2]; | |
760 | vV(0,2) = cov[3]; vV(2,0) = cov[3]; | |
761 | vV(1,2) = cov[4]; vV(2,1) = cov[4]; | |
762 | vV(2,2) = cov[5]; | |
763 | ||
764 | TMatrixD sumWi(TMatrixD::kInverted,vV); | |
765 | TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv); | |
766 | ||
dad616ce | 767 | Int_t nUsedTrks = inVtx->GetNIndices(); |
fee9ef0d | 768 | Double_t chi2 = inVtx->GetChi2(); |
769 | ||
6d8df534 | 770 | AliExternalTrackParam *track = 0; |
771 | Int_t ntrks = (Int_t)trkArray->GetEntriesFast(); | |
fee9ef0d | 772 | for(Int_t i=0;i<ntrks;i++) { |
6d8df534 | 773 | track = (AliExternalTrackParam*)trkArray->At(i); |
774 | if(!inVtx->UsesTrack(id[i])) { | |
919e537f | 775 | printf("track %d was not used in vertex fit",id[i]); |
fee9ef0d | 776 | continue; |
777 | } | |
778 | Double_t alpha = track->GetAlpha(); | |
779 | Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha); | |
6d8df534 | 780 | track->PropagateTo(xl,GetFieldkG()); |
fee9ef0d | 781 | // vector of track global coordinates |
782 | TMatrixD ri(3,1); | |
783 | // covariance matrix of ri | |
784 | TMatrixD wWi(3,3); | |
785 | ||
786 | // get space point from track | |
787 | if(!TrackToPoint(track,ri,wWi)) continue; | |
788 | ||
789 | TMatrixD wWiri(wWi,TMatrixD::kMult,ri); | |
790 | ||
791 | sumWi -= wWi; | |
792 | sumWiri -= wWiri; | |
793 | ||
6b908072 | 794 | // track contribution to chi2 |
fee9ef0d | 795 | TMatrixD deltar = rv; deltar -= ri; |
796 | TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar); | |
797 | Double_t chi2i = deltar(0,0)*wWideltar(0,0)+ | |
798 | deltar(1,0)*wWideltar(1,0)+ | |
799 | deltar(2,0)*wWideltar(2,0); | |
800 | // remove from total chi2 | |
801 | chi2 -= chi2i; | |
802 | ||
803 | nUsedTrks--; | |
804 | if(nUsedTrks<2) { | |
919e537f | 805 | printf("Trying to remove too many tracks!"); |
fee9ef0d | 806 | return 0x0; |
807 | } | |
808 | } | |
809 | ||
810 | TMatrixD rvnew(3,1); | |
811 | TMatrixD vVnew(3,3); | |
812 | ||
813 | // new inverted of weights matrix | |
814 | TMatrixD invsumWi(TMatrixD::kInverted,sumWi); | |
815 | vVnew = invsumWi; | |
816 | // new position of primary vertex | |
817 | rvnew.Mult(vVnew,sumWiri); | |
818 | ||
819 | Double_t position[3]; | |
820 | position[0] = rvnew(0,0); | |
821 | position[1] = rvnew(1,0); | |
822 | position[2] = rvnew(2,0); | |
823 | cov[0] = vVnew(0,0); | |
824 | cov[1] = vVnew(0,1); | |
825 | cov[2] = vVnew(1,1); | |
826 | cov[3] = vVnew(0,2); | |
827 | cov[4] = vVnew(1,2); | |
828 | cov[5] = vVnew(2,2); | |
829 | ||
830 | // store data in the vertex object | |
6b908072 | 831 | AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks+1); // the +1 is for the constraint |
fee9ef0d | 832 | outVtx->SetTitle(inVtx->GetTitle()); |
fee9ef0d | 833 | UShort_t *inindices = inVtx->GetIndices(); |
dad616ce | 834 | Int_t nIndices = nUsedTrks; |
4d023a8f | 835 | UShort_t *outindices = new UShort_t[nIndices]; |
fee9ef0d | 836 | Int_t j=0; |
fee9ef0d | 837 | for(Int_t k=0; k<inVtx->GetNIndices(); k++) { |
6d8df534 | 838 | Bool_t copyindex=kTRUE; |
fee9ef0d | 839 | for(Int_t l=0; l<ntrks; l++) { |
6d8df534 | 840 | if(inindices[k]==id[l]) copyindex=kFALSE; |
fee9ef0d | 841 | } |
842 | if(copyindex) { | |
843 | outindices[j] = inindices[k]; j++; | |
844 | } | |
845 | } | |
4d023a8f | 846 | outVtx->SetIndices(nIndices,outindices); |
4ce766eb | 847 | if (outindices) delete [] outindices; |
fee9ef0d | 848 | |
919e537f | 849 | /* |
850 | printf("Vertex before removing tracks:"); | |
fee9ef0d | 851 | inVtx->PrintStatus(); |
852 | inVtx->PrintIndices(); | |
919e537f | 853 | printf("Vertex after removing tracks:"); |
fee9ef0d | 854 | outVtx->PrintStatus(); |
855 | outVtx->PrintIndices(); | |
919e537f | 856 | */ |
fee9ef0d | 857 | |
858 | return outVtx; | |
859 | } | |
860 | //--------------------------------------------------------------------------- | |
6b908072 | 861 | AliESDVertex* AliVertexerTracks::RemoveConstraintFromVertex(AliESDVertex *inVtx, |
862 | Float_t *diamondxyz, | |
863 | Float_t *diamondcov) const | |
864 | { | |
865 | // | |
866 | // Removes diamond constraint from fit of inVtx | |
867 | // | |
868 | ||
869 | if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) { | |
870 | printf("ERROR: primary vertex has no constraint: cannot remove it\n"); | |
871 | return 0x0; | |
872 | } | |
873 | if(inVtx->GetNContributors()<3) { | |
874 | printf("ERROR: primary vertex has less than 2 tracks: cannot remove contraint\n"); | |
875 | return 0x0; | |
876 | } | |
877 | ||
878 | // diamond constraint | |
879 | TMatrixD vVb(3,3); | |
880 | vVb(0,0) = diamondcov[0]; | |
881 | vVb(0,1) = diamondcov[1]; | |
882 | vVb(0,2) = 0.; | |
883 | vVb(1,0) = diamondcov[1]; | |
884 | vVb(1,1) = diamondcov[2]; | |
885 | vVb(1,2) = 0.; | |
886 | vVb(2,0) = 0.; | |
887 | vVb(2,1) = 0.; | |
888 | vVb(2,2) = diamondcov[5]; | |
889 | TMatrixD vVbInv(TMatrixD::kInverted,vVb); | |
890 | TMatrixD rb(3,1); | |
891 | rb(0,0) = diamondxyz[0]; | |
892 | rb(1,0) = diamondxyz[1]; | |
893 | rb(2,0) = diamondxyz[2]; | |
894 | TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb); | |
895 | ||
896 | // input vertex | |
897 | TMatrixD rv(3,1); | |
898 | rv(0,0) = inVtx->GetXv(); | |
899 | rv(1,0) = inVtx->GetYv(); | |
900 | rv(2,0) = inVtx->GetZv(); | |
901 | TMatrixD vV(3,3); | |
902 | Double_t cov[6]; | |
903 | inVtx->GetCovMatrix(cov); | |
904 | vV(0,0) = cov[0]; | |
905 | vV(0,1) = cov[1]; vV(1,0) = cov[1]; | |
906 | vV(1,1) = cov[2]; | |
907 | vV(0,2) = cov[3]; vV(2,0) = cov[3]; | |
908 | vV(1,2) = cov[4]; vV(2,1) = cov[4]; | |
909 | vV(2,2) = cov[5]; | |
910 | TMatrixD vVInv(TMatrixD::kInverted,vV); | |
911 | TMatrixD vVInvrv(vVInv,TMatrixD::kMult,rv); | |
912 | ||
913 | ||
914 | TMatrixD sumWi = vVInv - vVbInv; | |
915 | ||
916 | ||
917 | TMatrixD sumWiri = vVInvrv - vVbInvrb; | |
918 | ||
919 | TMatrixD rvnew(3,1); | |
920 | TMatrixD vVnew(3,3); | |
921 | ||
922 | // new inverted of weights matrix | |
923 | TMatrixD invsumWi(TMatrixD::kInverted,sumWi); | |
924 | vVnew = invsumWi; | |
925 | // new position of primary vertex | |
926 | rvnew.Mult(vVnew,sumWiri); | |
927 | ||
928 | Double_t position[3]; | |
929 | position[0] = rvnew(0,0); | |
930 | position[1] = rvnew(1,0); | |
931 | position[2] = rvnew(2,0); | |
932 | cov[0] = vVnew(0,0); | |
933 | cov[1] = vVnew(0,1); | |
934 | cov[2] = vVnew(1,1); | |
935 | cov[3] = vVnew(0,2); | |
936 | cov[4] = vVnew(1,2); | |
937 | cov[5] = vVnew(2,2); | |
938 | ||
939 | ||
940 | Double_t chi2 = inVtx->GetChi2(); | |
941 | ||
942 | // diamond constribution to chi2 | |
943 | TMatrixD deltar = rv; deltar -= rb; | |
944 | TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar); | |
945 | Double_t chi2b = deltar(0,0)*vVbInvdeltar(0,0)+ | |
946 | deltar(1,0)*vVbInvdeltar(1,0)+ | |
947 | deltar(2,0)*vVbInvdeltar(2,0); | |
948 | // remove from total chi2 | |
949 | chi2 -= chi2b; | |
950 | ||
951 | // store data in the vertex object | |
952 | AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,inVtx->GetNContributors()-1); | |
953 | outVtx->SetTitle("VertexerTracksNoConstraint"); | |
954 | UShort_t *inindices = inVtx->GetIndices(); | |
955 | Int_t nIndices = inVtx->GetNIndices(); | |
956 | outVtx->SetIndices(nIndices,inindices); | |
957 | ||
958 | return outVtx; | |
959 | } | |
960 | //--------------------------------------------------------------------------- | |
a00021a7 | 961 | void AliVertexerTracks::SetCuts(Double_t *cuts) |
962 | { | |
963 | // | |
964 | // Cut values | |
965 | // | |
966 | SetDCAcut(cuts[0]); | |
967 | SetDCAcutIter0(cuts[1]); | |
968 | SetMaxd0z0(cuts[2]); | |
5b78b791 | 969 | if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired(); |
970 | SetMinClusters((Int_t)(TMath::Abs(cuts[3]))); | |
a00021a7 | 971 | SetMinTracks((Int_t)(cuts[4])); |
972 | SetNSigmad0(cuts[5]); | |
973 | SetMinDetFitter(cuts[6]); | |
974 | SetMaxTgl(cuts[7]); | |
975 | SetFiducialRZ(cuts[8],cuts[9]); | |
8c75f668 | 976 | fAlgo=(Int_t)(cuts[10]); |
977 | fAlgoIter0=(Int_t)(cuts[11]); | |
a00021a7 | 978 | |
979 | return; | |
980 | } | |
981 | //--------------------------------------------------------------------------- | |
f09c879d | 982 | void AliVertexerTracks::SetITSMode(Double_t dcacut, |
983 | Double_t dcacutIter0, | |
984 | Double_t maxd0z0, | |
6b29d399 | 985 | Int_t minCls, |
f09c879d | 986 | Int_t mintrks, |
987 | Double_t nsigma, | |
988 | Double_t mindetfitter, | |
6b29d399 | 989 | Double_t maxtgl, |
990 | Double_t fidR, | |
8c75f668 | 991 | Double_t fidZ, |
992 | Int_t finderAlgo, | |
993 | Int_t finderAlgoIter0) | |
f09c879d | 994 | { |
995 | // | |
996 | // Cut values for ITS mode | |
997 | // | |
998 | fMode = 0; | |
5b78b791 | 999 | if(minCls>0) { |
1000 | SetITSrefitRequired(); | |
1001 | } else { | |
1002 | SetITSrefitNotRequired(); | |
1003 | } | |
f09c879d | 1004 | SetDCAcut(dcacut); |
1005 | SetDCAcutIter0(dcacutIter0); | |
1006 | SetMaxd0z0(maxd0z0); | |
5b78b791 | 1007 | SetMinClusters(TMath::Abs(minCls)); |
f09c879d | 1008 | SetMinTracks(mintrks); |
1009 | SetNSigmad0(nsigma); | |
1010 | SetMinDetFitter(mindetfitter); | |
1011 | SetMaxTgl(maxtgl); | |
6b29d399 | 1012 | SetFiducialRZ(fidR,fidZ); |
8c75f668 | 1013 | fAlgo=finderAlgo; |
1014 | fAlgoIter0=finderAlgoIter0; | |
f09c879d | 1015 | |
1016 | return; | |
1017 | } | |
1018 | //--------------------------------------------------------------------------- | |
1019 | void AliVertexerTracks::SetTPCMode(Double_t dcacut, | |
1020 | Double_t dcacutIter0, | |
1021 | Double_t maxd0z0, | |
6b29d399 | 1022 | Int_t minCls, |
f09c879d | 1023 | Int_t mintrks, |
1024 | Double_t nsigma, | |
1025 | Double_t mindetfitter, | |
6b29d399 | 1026 | Double_t maxtgl, |
1027 | Double_t fidR, | |
8c75f668 | 1028 | Double_t fidZ, |
1029 | Int_t finderAlgo, | |
1030 | Int_t finderAlgoIter0) | |
f09c879d | 1031 | { |
1032 | // | |
1033 | // Cut values for TPC mode | |
1034 | // | |
1035 | fMode = 1; | |
1036 | SetITSrefitNotRequired(); | |
1037 | SetDCAcut(dcacut); | |
1038 | SetDCAcutIter0(dcacutIter0); | |
1039 | SetMaxd0z0(maxd0z0); | |
6b29d399 | 1040 | SetMinClusters(minCls); |
f09c879d | 1041 | SetMinTracks(mintrks); |
1042 | SetNSigmad0(nsigma); | |
1043 | SetMinDetFitter(mindetfitter); | |
1044 | SetMaxTgl(maxtgl); | |
6b29d399 | 1045 | SetFiducialRZ(fidR,fidZ); |
8c75f668 | 1046 | fAlgo=finderAlgo; |
1047 | fAlgoIter0=finderAlgoIter0; | |
f09c879d | 1048 | |
1049 | return; | |
1050 | } | |
1051 | //--------------------------------------------------------------------------- | |
f5740e9a | 1052 | void AliVertexerTracks::SetSkipTracks(Int_t n,const Int_t *skipped) |
fee9ef0d | 1053 | { |
1054 | // | |
1055 | // Mark the tracks not to be used in the vertex reconstruction. | |
1056 | // Tracks are identified by AliESDtrack::GetID() | |
1057 | // | |
1058 | fNTrksToSkip = n; fTrksToSkip = new Int_t[n]; | |
50ff0bcd | 1059 | for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; |
1060 | return; | |
1061 | } | |
1062 | //--------------------------------------------------------------------------- | |
fee9ef0d | 1063 | void AliVertexerTracks::SetVtxStart(AliESDVertex *vtx) |
1064 | { | |
50ff0bcd | 1065 | // |
1066 | // Set initial vertex knowledge | |
1067 | // | |
1068 | vtx->GetXYZ(fNominalPos); | |
1069 | vtx->GetCovMatrix(fNominalCov); | |
fee9ef0d | 1070 | SetConstraintOn(); |
50ff0bcd | 1071 | return; |
1072 | } | |
2d57349e | 1073 | //--------------------------------------------------------------------------- |
fee9ef0d | 1074 | void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights) |
dc3719f3 | 1075 | { |
6d8df534 | 1076 | AliExternalTrackParam *track1; |
6d8df534 | 1077 | const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast(); |
8e104736 | 1078 | AliStrLine **linarray = new AliStrLine* [knacc]; |
dc3719f3 | 1079 | for(Int_t i=0; i<knacc; i++){ |
6d8df534 | 1080 | track1 = (AliExternalTrackParam*)fTrkArraySel.At(i); |
dc3719f3 | 1081 | Double_t alpha=track1->GetAlpha(); |
1082 | Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1]; | |
1083 | Double_t pos[3],dir[3],sigmasq[3]; | |
f09c879d | 1084 | track1->GetXYZAt(mindist,GetFieldkG(),pos); |
1085 | track1->GetPxPyPzAt(mindist,GetFieldkG(),dir); | |
dc3719f3 | 1086 | sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2(); |
1087 | sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2(); | |
1088 | sigmasq[2]=track1->GetSigmaZ2(); | |
146c29df | 1089 | TMatrixD ri(3,1); |
1090 | TMatrixD wWi(3,3); | |
8c75f668 | 1091 | if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");} |
146c29df | 1092 | Double_t wmat[9]; |
1093 | Int_t iel=0; | |
1094 | for(Int_t ia=0;ia<3;ia++){ | |
1095 | for(Int_t ib=0;ib<3;ib++){ | |
1096 | wmat[iel]=wWi(ia,ib); | |
1097 | iel++; | |
1098 | } | |
1099 | } | |
8e104736 | 1100 | linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir); |
dc3719f3 | 1101 | } |
8e104736 | 1102 | fVert=TrackletVertexFinder(linarray,knacc,optUseWeights); |
1103 | for(Int_t i=0; i<knacc; i++) delete linarray[i]; | |
1104 | delete [] linarray; | |
dc3719f3 | 1105 | } |
1106 | //--------------------------------------------------------------------------- | |
f5740e9a | 1107 | AliESDVertex AliVertexerTracks::TrackletVertexFinder(const TClonesArray *lines, Int_t optUseWeights) |
fee9ef0d | 1108 | { |
8e104736 | 1109 | // Calculate the point at minimum distance to prepared tracks (TClonesArray) |
dc3719f3 | 1110 | const Int_t knacc = (Int_t)lines->GetEntriesFast(); |
8e104736 | 1111 | AliStrLine** lines2 = new AliStrLine* [knacc]; |
1112 | for(Int_t i=0; i<knacc; i++){ | |
1113 | lines2[i]= (AliStrLine*)lines->At(i); | |
1114 | } | |
1115 | AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights); | |
1116 | delete [] lines2; | |
1117 | return vert; | |
1118 | } | |
1119 | ||
1120 | //--------------------------------------------------------------------------- | |
1121 | AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights) | |
1122 | { | |
1123 | // Calculate the point at minimum distance to prepared tracks (array of AliStrLine) | |
1124 | ||
dc3719f3 | 1125 | Double_t initPos[3]={0.,0.,0.}; |
2d57349e | 1126 | |
9d345960 | 1127 | Double_t (*vectP0)[3]=new Double_t [knacc][3](); |
1128 | Double_t (*vectP1)[3]=new Double_t [knacc][3](); | |
2d57349e | 1129 | |
1130 | Double_t sum[3][3]; | |
1131 | Double_t dsum[3]={0,0,0}; | |
146c29df | 1132 | TMatrixD sumWi(3,3); |
1133 | for(Int_t i=0;i<3;i++){ | |
1134 | for(Int_t j=0;j<3;j++){ | |
1135 | sum[i][j]=0; | |
1136 | sumWi(i,j)=0.; | |
1137 | } | |
1138 | } | |
1139 | ||
2d57349e | 1140 | for(Int_t i=0; i<knacc; i++){ |
8e104736 | 1141 | AliStrLine *line1 = lines[i]; |
dc3719f3 | 1142 | Double_t p0[3],cd[3],sigmasq[3]; |
146c29df | 1143 | Double_t wmat[9]; |
27275ff3 | 1144 | if(!line1) { printf("ERROR %d %d\n",i,knacc); continue; } |
2d57349e | 1145 | line1->GetP0(p0); |
1146 | line1->GetCd(cd); | |
dc3719f3 | 1147 | line1->GetSigma2P0(sigmasq); |
146c29df | 1148 | line1->GetWMatrix(wmat); |
1149 | TMatrixD wWi(3,3); | |
1150 | Int_t iel=0; | |
1151 | for(Int_t ia=0;ia<3;ia++){ | |
1152 | for(Int_t ib=0;ib<3;ib++){ | |
1153 | wWi(ia,ib)=wmat[iel]; | |
1154 | iel++; | |
1155 | } | |
1156 | } | |
1157 | ||
1158 | sumWi+=wWi; | |
1159 | ||
2d57349e | 1160 | Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]}; |
1161 | vectP0[i][0]=p0[0]; | |
1162 | vectP0[i][1]=p0[1]; | |
1163 | vectP0[i][2]=p0[2]; | |
1164 | vectP1[i][0]=p1[0]; | |
1165 | vectP1[i][1]=p1[1]; | |
1166 | vectP1[i][2]=p1[2]; | |
1167 | ||
1168 | Double_t matr[3][3]; | |
1169 | Double_t dknow[3]; | |
1170 | if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow); | |
dc3719f3 | 1171 | else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow); |
1172 | ||
2d57349e | 1173 | |
1174 | for(Int_t iii=0;iii<3;iii++){ | |
1175 | dsum[iii]+=dknow[iii]; | |
1176 | for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj]; | |
1177 | } | |
2d57349e | 1178 | } |
146c29df | 1179 | |
1180 | TMatrixD invsumWi(TMatrixD::kInverted,sumWi); | |
1181 | Double_t covmatrix[6]; | |
1182 | covmatrix[0] = invsumWi(0,0); | |
1183 | covmatrix[1] = invsumWi(0,1); | |
1184 | covmatrix[2] = invsumWi(1,1); | |
1185 | covmatrix[3] = invsumWi(0,2); | |
1186 | covmatrix[4] = invsumWi(1,2); | |
1187 | covmatrix[5] = invsumWi(2,2); | |
1188 | ||
2d57349e | 1189 | Double_t vett[3][3]; |
1190 | Double_t det=GetDeterminant3X3(sum); | |
dc3719f3 | 1191 | Double_t sigma=0; |
2d57349e | 1192 | |
f5740e9a | 1193 | if(TMath::Abs(det) > kAlmost0){ |
fee9ef0d | 1194 | for(Int_t zz=0;zz<3;zz++){ |
1195 | for(Int_t ww=0;ww<3;ww++){ | |
1196 | for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk]; | |
1197 | } | |
1198 | for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk]; | |
1199 | initPos[zz]=GetDeterminant3X3(vett)/det; | |
1200 | } | |
1201 | ||
1202 | ||
1203 | for(Int_t i=0; i<knacc; i++){ | |
1204 | Double_t p0[3]={0,0,0},p1[3]={0,0,0}; | |
1205 | for(Int_t ii=0;ii<3;ii++){ | |
1206 | p0[ii]=vectP0[i][ii]; | |
1207 | p1[ii]=vectP1[i][ii]; | |
1208 | } | |
1209 | sigma+=GetStrLinMinDist(p0,p1,initPos); | |
1210 | } | |
1211 | ||
f09c879d | 1212 | if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;} |
fee9ef0d | 1213 | }else{ |
50ff0bcd | 1214 | sigma=999; |
1215 | } | |
146c29df | 1216 | AliESDVertex theVert(initPos,covmatrix,99999.,knacc); |
1217 | theVert.SetDispersion(sigma); | |
4ce766eb | 1218 | delete [] vectP0; |
1219 | delete [] vectP1; | |
dc3719f3 | 1220 | return theVert; |
50ff0bcd | 1221 | } |
8e104736 | 1222 | |
50ff0bcd | 1223 | //--------------------------------------------------------------------------- |
6d8df534 | 1224 | Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t, |
817e1004 | 1225 | TMatrixD &ri,TMatrixD &wWi, |
0bf2e2b4 | 1226 | Bool_t uUi3by3) const |
fee9ef0d | 1227 | { |
1228 | // | |
6d8df534 | 1229 | // Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix |
fee9ef0d | 1230 | // wWi of the space point that it represents (to be used in VertexFitter()) |
1231 | // | |
1232 | ||
1233 | ||
1234 | Double_t rotAngle = t->GetAlpha(); | |
1235 | if(rotAngle<0.) rotAngle += 2.*TMath::Pi(); | |
1236 | Double_t cosRot = TMath::Cos(rotAngle); | |
1237 | Double_t sinRot = TMath::Sin(rotAngle); | |
1238 | ||
1239 | ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot; | |
1240 | ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot; | |
1241 | ri(2,0) = t->GetZ(); | |
1242 | ||
0bf2e2b4 | 1243 | |
1244 | ||
1245 | if(!uUi3by3) { | |
1246 | // matrix to go from global (x,y,z) to local (y,z); | |
1247 | TMatrixD qQi(2,3); | |
1248 | qQi(0,0) = -sinRot; | |
1249 | qQi(0,1) = cosRot; | |
1250 | qQi(0,2) = 0.; | |
1251 | qQi(1,0) = 0.; | |
1252 | qQi(1,1) = 0.; | |
1253 | qQi(1,2) = 1.; | |
1254 | ||
1255 | // covariance matrix of local (y,z) - inverted | |
0bf2e2b4 | 1256 | TMatrixD uUi(2,2); |
6d8df534 | 1257 | uUi(0,0) = t->GetSigmaY2(); |
1258 | uUi(0,1) = t->GetSigmaZY(); | |
1259 | uUi(1,0) = t->GetSigmaZY(); | |
1260 | uUi(1,1) = t->GetSigmaZ2(); | |
919e537f | 1261 | //printf(" Ui :"); |
1262 | //printf(" %f %f",uUi(0,0),uUi(0,1)); | |
1263 | //printf(" %f %f",uUi(1,0),uUi(1,1)); | |
0bf2e2b4 | 1264 | |
1265 | if(uUi.Determinant() <= 0.) return kFALSE; | |
1266 | TMatrixD uUiInv(TMatrixD::kInverted,uUi); | |
1267 | ||
1268 | // weights matrix: wWi = qQiT * uUiInv * qQi | |
1269 | TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi); | |
1270 | TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi); | |
1271 | wWi = m; | |
1272 | } else { | |
1273 | if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty"); | |
1274 | // matrix to go from global (x,y,z) to local (x,y,z); | |
1275 | TMatrixD qQi(3,3); | |
1276 | qQi(0,0) = cosRot; | |
1277 | qQi(0,1) = sinRot; | |
1278 | qQi(0,2) = 0.; | |
1279 | qQi(1,0) = -sinRot; | |
1280 | qQi(1,1) = cosRot; | |
1281 | qQi(1,2) = 0.; | |
1282 | qQi(2,0) = 0.; | |
1283 | qQi(2,1) = 0.; | |
1284 | qQi(2,2) = 1.; | |
1285 | ||
1286 | // covariance of fVert along the track | |
1287 | Double_t p[3],pt,ptot; | |
1288 | t->GetPxPyPz(p); | |
1289 | pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]); | |
1290 | ptot = TMath::Sqrt(pt*pt+p[2]*p[2]); | |
3b768a1d | 1291 | Double_t cphi = p[0]/pt; //cos(phi)=px/pt |
1292 | Double_t sphi = p[1]/pt; //sin(phi)=py/pt | |
1293 | Double_t clambda = pt/ptot; //cos(lambda)=pt/ptot | |
1294 | Double_t slambda = p[2]/ptot; //sin(lambda)=pz/ptot | |
0bf2e2b4 | 1295 | Double_t covfVert[6]; |
1296 | fVert.GetCovMatrix(covfVert); | |
1297 | Double_t covfVertalongt = | |
3b768a1d | 1298 | covfVert[0]*cphi*cphi*clambda*clambda |
1299 | +covfVert[1]*2.*cphi*sphi*clambda*clambda | |
1300 | +covfVert[3]*2.*cphi*clambda*slambda | |
1301 | +covfVert[2]*sphi*sphi*clambda*clambda | |
1302 | +covfVert[4]*2.*sphi*clambda*slambda | |
1303 | +covfVert[5]*slambda*slambda; | |
0bf2e2b4 | 1304 | // covariance matrix of local (x,y,z) - inverted |
1305 | TMatrixD uUi(3,3); | |
3b768a1d | 1306 | uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00; |
919e537f | 1307 | AliDebug(1,Form("=====> sqrtUi00 cm %f",TMath::Sqrt(uUi(0,0)))); |
0bf2e2b4 | 1308 | uUi(0,1) = 0.; |
1309 | uUi(0,2) = 0.; | |
1310 | uUi(1,0) = 0.; | |
6d8df534 | 1311 | uUi(1,1) = t->GetSigmaY2(); |
1312 | uUi(1,2) = t->GetSigmaZY(); | |
0bf2e2b4 | 1313 | uUi(2,0) = 0.; |
6d8df534 | 1314 | uUi(2,1) = t->GetSigmaZY(); |
1315 | uUi(2,2) = t->GetSigmaZ2(); | |
0bf2e2b4 | 1316 | //printf(" Ui :\n"); |
1317 | //printf(" %f %f\n",uUi(0,0),uUi(0,1)); | |
1318 | //printf(" %f %f\n",uUi(1,0),uUi(1,1)); | |
fee9ef0d | 1319 | |
0bf2e2b4 | 1320 | if(uUi.Determinant() <= 0.) return kFALSE; |
1321 | TMatrixD uUiInv(TMatrixD::kInverted,uUi); | |
1322 | ||
1323 | // weights matrix: wWi = qQiT * uUiInv * qQi | |
1324 | TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi); | |
1325 | TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi); | |
1326 | wWi = m; | |
1327 | } | |
1328 | ||
fee9ef0d | 1329 | |
1330 | return kTRUE; | |
1331 | } | |
1332 | //--------------------------------------------------------------------------- | |
6d8df534 | 1333 | void AliVertexerTracks::TooFewTracks() |
fee9ef0d | 1334 | { |
50ff0bcd | 1335 | // |
6d8df534 | 1336 | // When the number of tracks is < fMinTracks, |
1337 | // deal with vertices not found and prepare to exit | |
50ff0bcd | 1338 | // |
919e537f | 1339 | AliDebug(1,"TooFewTracks"); |
2d57349e | 1340 | |
50ff0bcd | 1341 | Double_t pos[3],err[3]; |
50ff0bcd | 1342 | pos[0] = fNominalPos[0]; |
1343 | err[0] = TMath::Sqrt(fNominalCov[0]); | |
1344 | pos[1] = fNominalPos[1]; | |
1345 | err[1] = TMath::Sqrt(fNominalCov[2]); | |
6d8df534 | 1346 | pos[2] = fNominalPos[2]; |
1347 | err[2] = TMath::Sqrt(fNominalCov[5]); | |
1348 | Int_t ncontr = (err[0]>1. ? -1 : -3); | |
5b78b791 | 1349 | if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; } |
50ff0bcd | 1350 | fCurrentVertex = new AliESDVertex(pos,err); |
1351 | fCurrentVertex->SetNContributors(ncontr); | |
2d57349e | 1352 | |
fee9ef0d | 1353 | if(fConstraint) { |
1354 | fCurrentVertex->SetTitle("VertexerTracksWithConstraint"); | |
1355 | } else { | |
50ff0bcd | 1356 | fCurrentVertex->SetTitle("VertexerTracksNoConstraint"); |
fee9ef0d | 1357 | } |
2d57349e | 1358 | |
6d8df534 | 1359 | if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete(); |
6fefa5ce | 1360 | if(fIdSel) {delete [] fIdSel; fIdSel=NULL;} |
1361 | if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;} | |
6d8df534 | 1362 | |
50ff0bcd | 1363 | return; |
1364 | } | |
1365 | //--------------------------------------------------------------------------- | |
fee9ef0d | 1366 | void AliVertexerTracks::VertexFinder(Int_t optUseWeights) |
1367 | { | |
2d57349e | 1368 | |
50ff0bcd | 1369 | // Get estimate of vertex position in (x,y) from tracks DCA |
1370 | ||
1371 | Double_t initPos[3]; | |
1372 | initPos[2] = 0.; | |
1373 | for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i]; | |
6d8df534 | 1374 | Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast(); |
50ff0bcd | 1375 | Double_t aver[3]={0.,0.,0.}; |
1376 | Double_t aversq[3]={0.,0.,0.}; | |
1377 | Double_t sigmasq[3]={0.,0.,0.}; | |
1378 | Double_t sigma=0; | |
1379 | Int_t ncombi = 0; | |
6d8df534 | 1380 | AliExternalTrackParam *track1; |
1381 | AliExternalTrackParam *track2; | |
50ff0bcd | 1382 | Double_t pos[3],dir[3]; |
1383 | Double_t alpha,mindist; | |
2d57349e | 1384 | |
50ff0bcd | 1385 | for(Int_t i=0; i<nacc; i++){ |
6d8df534 | 1386 | track1 = (AliExternalTrackParam*)fTrkArraySel.At(i); |
50ff0bcd | 1387 | alpha=track1->GetAlpha(); |
1388 | mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1]; | |
f09c879d | 1389 | track1->GetXYZAt(mindist,GetFieldkG(),pos); |
1390 | track1->GetPxPyPzAt(mindist,GetFieldkG(),dir); | |
50ff0bcd | 1391 | AliStrLine *line1 = new AliStrLine(pos,dir); |
2d57349e | 1392 | |
50ff0bcd | 1393 | // AliStrLine *line1 = new AliStrLine(); |
f09c879d | 1394 | // track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1); |
50ff0bcd | 1395 | |
1396 | for(Int_t j=i+1; j<nacc; j++){ | |
6d8df534 | 1397 | track2 = (AliExternalTrackParam*)fTrkArraySel.At(j); |
50ff0bcd | 1398 | alpha=track2->GetAlpha(); |
1399 | mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1]; | |
f09c879d | 1400 | track2->GetXYZAt(mindist,GetFieldkG(),pos); |
1401 | track2->GetPxPyPzAt(mindist,GetFieldkG(),dir); | |
50ff0bcd | 1402 | AliStrLine *line2 = new AliStrLine(pos,dir); |
1403 | // AliStrLine *line2 = new AliStrLine(); | |
f09c879d | 1404 | // track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2); |
50ff0bcd | 1405 | Double_t distCA=line2->GetDCA(line1); |
fee9ef0d | 1406 | //printf("%d %d %f\n",i,j,distCA); |
1407 | if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){ | |
50ff0bcd | 1408 | Double_t pnt1[3],pnt2[3],crosspoint[3]; |
2d57349e | 1409 | |
50ff0bcd | 1410 | if(optUseWeights<=0){ |
1411 | Int_t retcode = line2->Cross(line1,crosspoint); | |
1412 | if(retcode>=0){ | |
1413 | ncombi++; | |
1414 | for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj]; | |
1415 | for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]); | |
1416 | } | |
1417 | } | |
1418 | if(optUseWeights>0){ | |
1419 | Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2); | |
1420 | if(retcode>=0){ | |
5c03260f | 1421 | Double_t cs, sn; |
50ff0bcd | 1422 | alpha=track1->GetAlpha(); |
1423 | cs=TMath::Cos(alpha); sn=TMath::Sin(alpha); | |
1424 | Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2(); | |
1425 | alpha=track2->GetAlpha(); | |
1426 | cs=TMath::Cos(alpha); sn=TMath::Sin(alpha); | |
1427 | Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2(); | |
1428 | Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2(); | |
1429 | Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1; | |
1430 | Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1; | |
1431 | Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1; | |
1432 | crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0]; | |
1433 | crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1]; | |
1434 | crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2]; | |
1435 | ||
1436 | ncombi++; | |
1437 | for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj]; | |
1438 | for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]); | |
1439 | } | |
1440 | } | |
1441 | } | |
1442 | delete line2; | |
1443 | } | |
1444 | delete line1; | |
71d84967 | 1445 | } |
50ff0bcd | 1446 | if(ncombi>0){ |
1447 | for(Int_t jj=0;jj<3;jj++){ | |
1448 | initPos[jj] = aver[jj]/ncombi; | |
fee9ef0d | 1449 | //printf("%f\n",initPos[jj]); |
50ff0bcd | 1450 | aversq[jj]/=ncombi; |
1451 | sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj]; | |
1452 | sigma+=sigmasq[jj]; | |
1453 | } | |
1454 | sigma=TMath::Sqrt(TMath::Abs(sigma)); | |
1455 | } | |
1456 | else { | |
1457 | Warning("VertexFinder","Finder did not succed"); | |
1458 | sigma=999; | |
1459 | } | |
1460 | fVert.SetXYZ(initPos); | |
1461 | fVert.SetDispersion(sigma); | |
1462 | fVert.SetNContributors(ncombi); | |
2d57349e | 1463 | } |
ec1be5d5 | 1464 | //--------------------------------------------------------------------------- |
6d8df534 | 1465 | void AliVertexerTracks::VertexFitter() |
fee9ef0d | 1466 | { |
ec1be5d5 | 1467 | // |
1468 | // The optimal estimate of the vertex position is given by a "weighted | |
0bf2e2b4 | 1469 | // average of tracks positions". |
fee9ef0d | 1470 | // Original method: V. Karimaki, CMS Note 97/0051 |
ec1be5d5 | 1471 | // |
6d8df534 | 1472 | Bool_t useConstraint = fConstraint; |
ec1be5d5 | 1473 | Double_t initPos[3]; |
0bf2e2b4 | 1474 | if(!fOnlyFitter) { |
1475 | fVert.GetXYZ(initPos); | |
1476 | } else { | |
1477 | initPos[0]=fNominalPos[0]; | |
1478 | initPos[1]=fNominalPos[1]; | |
1479 | initPos[2]=fNominalPos[2]; | |
1480 | } | |
1481 | ||
6d8df534 | 1482 | Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries(); |
1483 | if(nTrksSel==1) useConstraint=kTRUE; | |
919e537f | 1484 | AliDebug(1,"--- VertexFitter(): start"); |
1485 | AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel)); | |
1486 | AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks)); | |
1487 | AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2])); | |
1488 | if(useConstraint) AliDebug(1,Form(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],TMath::Sqrt(fNominalCov[0]),fNominalPos[1],TMath::Sqrt(fNominalCov[2]))); | |
ec1be5d5 | 1489 | |
0bf2e2b4 | 1490 | // special treatment for few-tracks fits (e.g. secondary vertices) |
6d8df534 | 1491 | Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE; |
0bf2e2b4 | 1492 | |
ec1be5d5 | 1493 | Int_t i,j,k,step=0; |
1494 | TMatrixD rv(3,1); | |
1495 | TMatrixD vV(3,3); | |
1496 | rv(0,0) = initPos[0]; | |
1497 | rv(1,0) = initPos[1]; | |
1498 | rv(2,0) = 0.; | |
1499 | Double_t xlStart,alpha; | |
6d8df534 | 1500 | Int_t nTrksUsed; |
fee9ef0d | 1501 | Double_t chi2,chi2i,chi2b; |
6d8df534 | 1502 | AliExternalTrackParam *t = 0; |
ec1be5d5 | 1503 | Int_t failed = 0; |
1504 | ||
50ff0bcd | 1505 | // initial vertex covariance matrix |
1506 | TMatrixD vVb(3,3); | |
1507 | vVb(0,0) = fNominalCov[0]; | |
1508 | vVb(0,1) = fNominalCov[1]; | |
1509 | vVb(0,2) = 0.; | |
1510 | vVb(1,0) = fNominalCov[1]; | |
1511 | vVb(1,1) = fNominalCov[2]; | |
1512 | vVb(1,2) = 0.; | |
1513 | vVb(2,0) = 0.; | |
1514 | vVb(2,1) = 0.; | |
1515 | vVb(2,2) = fNominalCov[5]; | |
1516 | TMatrixD vVbInv(TMatrixD::kInverted,vVb); | |
1517 | TMatrixD rb(3,1); | |
1518 | rb(0,0) = fNominalPos[0]; | |
1519 | rb(1,0) = fNominalPos[1]; | |
1520 | rb(2,0) = fNominalPos[2]; | |
1521 | TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb); | |
1522 | ||
1523 | ||
eae1677e | 1524 | // 2 steps: |
1525 | // 1st - estimate of vtx using all tracks | |
1526 | // 2nd - estimate of global chi2 | |
1527 | for(step=0; step<2; step++) { | |
919e537f | 1528 | AliDebug(1,Form(" step = %d\n",step)); |
ec1be5d5 | 1529 | chi2 = 0.; |
6d8df534 | 1530 | nTrksUsed = 0; |
ec1be5d5 | 1531 | |
0bf2e2b4 | 1532 | if(step==1) { initPos[0]=rv(0,0); initPos[0]=rv(1,0); } |
1533 | ||
ec1be5d5 | 1534 | TMatrixD sumWiri(3,1); |
1535 | TMatrixD sumWi(3,3); | |
1536 | for(i=0; i<3; i++) { | |
1537 | sumWiri(i,0) = 0.; | |
1538 | for(j=0; j<3; j++) sumWi(j,i) = 0.; | |
1539 | } | |
1540 | ||
fee9ef0d | 1541 | // mean vertex constraint |
1542 | if(useConstraint) { | |
50ff0bcd | 1543 | for(i=0;i<3;i++) { |
1544 | sumWiri(i,0) += vVbInvrb(i,0); | |
1545 | for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k); | |
07680cae | 1546 | } |
fee9ef0d | 1547 | // chi2 |
1548 | TMatrixD deltar = rv; deltar -= rb; | |
1549 | TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar); | |
1550 | chi2b = deltar(0,0)*vVbInvdeltar(0,0)+ | |
1551 | deltar(1,0)*vVbInvdeltar(1,0)+ | |
1552 | deltar(2,0)*vVbInvdeltar(2,0); | |
1553 | chi2 += chi2b; | |
07680cae | 1554 | } |
1555 | ||
ec1be5d5 | 1556 | // loop on tracks |
6d8df534 | 1557 | for(k=0; k<nTrksSel; k++) { |
ec1be5d5 | 1558 | // get track from track array |
6d8df534 | 1559 | t = (AliExternalTrackParam*)fTrkArraySel.At(k); |
ec1be5d5 | 1560 | alpha = t->GetAlpha(); |
1561 | xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha); | |
fee9ef0d | 1562 | // to vtxSeed (from finder) |
6d8df534 | 1563 | t->PropagateTo(xlStart,GetFieldkG()); |
fee9ef0d | 1564 | |
ec1be5d5 | 1565 | // vector of track global coordinates |
1566 | TMatrixD ri(3,1); | |
fee9ef0d | 1567 | // covariance matrix of ri |
1568 | TMatrixD wWi(3,3); | |
1569 | ||
1570 | // get space point from track | |
0bf2e2b4 | 1571 | if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue; |
fee9ef0d | 1572 | TMatrixD wWiri(wWi,TMatrixD::kMult,ri); |
ec1be5d5 | 1573 | |
1574 | // track chi2 | |
1575 | TMatrixD deltar = rv; deltar -= ri; | |
1576 | TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar); | |
1577 | chi2i = deltar(0,0)*wWideltar(0,0)+ | |
1578 | deltar(1,0)*wWideltar(1,0)+ | |
1579 | deltar(2,0)*wWideltar(2,0); | |
1580 | ||
ec1be5d5 | 1581 | // add to total chi2 |
1582 | chi2 += chi2i; | |
1583 | ||
ec1be5d5 | 1584 | sumWiri += wWiri; |
1585 | sumWi += wWi; | |
1586 | ||
6d8df534 | 1587 | nTrksUsed++; |
ec1be5d5 | 1588 | } // end loop on tracks |
1589 | ||
6d8df534 | 1590 | if(nTrksUsed < fMinTracks) { |
ec1be5d5 | 1591 | failed=1; |
1592 | continue; | |
1593 | } | |
1594 | ||
1595 | Double_t determinant = sumWi.Determinant(); | |
f09c879d | 1596 | if(determinant < fMinDetFitter) { |
919e537f | 1597 | AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter)); |
ec1be5d5 | 1598 | failed=1; |
1599 | continue; | |
1600 | } | |
1601 | ||
0bf2e2b4 | 1602 | if(step==0) { |
1603 | // inverted of weights matrix | |
1604 | TMatrixD invsumWi(TMatrixD::kInverted,sumWi); | |
1605 | vV = invsumWi; | |
1606 | // position of primary vertex | |
1607 | rv.Mult(vV,sumWiri); | |
1608 | } | |
eae1677e | 1609 | } // end loop on the 2 steps |
ec1be5d5 | 1610 | |
ec1be5d5 | 1611 | if(failed) { |
6d8df534 | 1612 | TooFewTracks(); |
ec1be5d5 | 1613 | return; |
1614 | } | |
1615 | ||
1616 | Double_t position[3]; | |
1617 | position[0] = rv(0,0); | |
1618 | position[1] = rv(1,0); | |
1619 | position[2] = rv(2,0); | |
1620 | Double_t covmatrix[6]; | |
1621 | covmatrix[0] = vV(0,0); | |
1622 | covmatrix[1] = vV(0,1); | |
1623 | covmatrix[2] = vV(1,1); | |
1624 | covmatrix[3] = vV(0,2); | |
1625 | covmatrix[4] = vV(1,2); | |
1626 | covmatrix[5] = vV(2,2); | |
1627 | ||
dc3719f3 | 1628 | // for correct chi2/ndf, count constraint as additional "track" |
6d8df534 | 1629 | if(fConstraint) nTrksUsed++; |
ec1be5d5 | 1630 | // store data in the vertex object |
5b78b791 | 1631 | if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; } |
6d8df534 | 1632 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed); |
ec1be5d5 | 1633 | |
919e537f | 1634 | AliDebug(1," Vertex after fit:"); |
1635 | AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors())); | |
1636 | AliDebug(1,"--- VertexFitter(): finish\n"); | |
1637 | ||
ec1be5d5 | 1638 | |
1639 | return; | |
1640 | } | |
50ff0bcd | 1641 | //---------------------------------------------------------------------------- |
f5740e9a | 1642 | AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArray, |
6d8df534 | 1643 | UShort_t *id, |
1644 | Bool_t optUseFitter, | |
97669588 | 1645 | Bool_t optPropagate, |
1646 | Bool_t optUseDiamondConstraint) | |
fee9ef0d | 1647 | { |
eae1677e | 1648 | // |
6d8df534 | 1649 | // Return vertex from tracks (AliExternalTrackParam) in array |
eae1677e | 1650 | // |
5b78b791 | 1651 | fCurrentVertex = 0; |
97669588 | 1652 | // set optUseDiamondConstraint=TRUE only if you are reconstructing the |
1653 | // primary vertex! | |
1654 | if(optUseDiamondConstraint) { | |
1655 | SetConstraintOn(); | |
1656 | } else { | |
1657 | SetConstraintOff(); | |
1658 | } | |
fee9ef0d | 1659 | |
50ff0bcd | 1660 | // get tracks and propagate them to initial vertex position |
6d8df534 | 1661 | fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()]; |
1662 | Int_t nTrksSel = PrepareTracks(*trkArray,id,0); | |
97669588 | 1663 | if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) || |
1664 | (optUseDiamondConstraint && nTrksSel<1)) { | |
6d8df534 | 1665 | TooFewTracks(); |
fee9ef0d | 1666 | return fCurrentVertex; |
50ff0bcd | 1667 | } |
1668 | ||
6d8df534 | 1669 | // vertex finder |
97669588 | 1670 | if(nTrksSel==1) { |
1671 | AliDebug(1,"Just one track"); | |
1672 | OneTrackVertFinder(); | |
1673 | } else { | |
1674 | switch (fAlgo) { | |
1675 | case 1: StrLinVertexFinderMinDist(1); break; | |
1676 | case 2: StrLinVertexFinderMinDist(0); break; | |
1677 | case 3: HelixVertexFinder(); break; | |
1678 | case 4: VertexFinder(1); break; | |
1679 | case 5: VertexFinder(0); break; | |
1680 | default: printf("Wrong algorithm\n"); break; | |
1681 | } | |
fee9ef0d | 1682 | } |
919e537f | 1683 | AliDebug(1," Vertex finding completed\n"); |
fee9ef0d | 1684 | |
1685 | // vertex fitter | |
6d8df534 | 1686 | if(optUseFitter) { |
1687 | VertexFitter(); | |
1688 | } else { | |
fee9ef0d | 1689 | Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()}; |
146c29df | 1690 | Double_t covmatrix[6]; |
1691 | fVert.GetCovMatrix(covmatrix); | |
fee9ef0d | 1692 | Double_t chi2=99999.; |
6d8df534 | 1693 | Int_t nTrksUsed=fVert.GetNContributors(); |
1694 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed); | |
fee9ef0d | 1695 | } |
1696 | fCurrentVertex->SetDispersion(fVert.GetDispersion()); | |
1697 | ||
1698 | ||
1699 | // set indices of used tracks and propagate track to found vertex | |
50ff0bcd | 1700 | UShort_t *indices = 0; |
6d8df534 | 1701 | Double_t d0z0[2],covd0z0[3]; |
1702 | AliExternalTrackParam *t = 0; | |
fee9ef0d | 1703 | if(fCurrentVertex->GetNContributors()>0) { |
44a76311 | 1704 | indices = new UShort_t[fTrkArraySel.GetEntriesFast()]; |
6d8df534 | 1705 | for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) { |
1706 | indices[jj] = fIdSel[jj]; | |
1707 | t = (AliExternalTrackParam*)fTrkArraySel.At(jj); | |
1708 | if(optPropagate && optUseFitter) { | |
fee9ef0d | 1709 | if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) { |
6d8df534 | 1710 | t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0); |
919e537f | 1711 | AliDebug(1,Form("Track %d propagated to found vertex",jj)); |
6d8df534 | 1712 | } else { |
fee9ef0d | 1713 | AliWarning("Found vertex outside beam pipe!"); |
1714 | } | |
1715 | } | |
50ff0bcd | 1716 | } |
fee9ef0d | 1717 | fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices); |
50ff0bcd | 1718 | } |
6d8df534 | 1719 | |
1720 | // clean up | |
4ce766eb | 1721 | if (indices) {delete [] indices; indices=NULL;} |
6fefa5ce | 1722 | delete [] fIdSel; fIdSel=NULL; |
6d8df534 | 1723 | fTrkArraySel.Delete(); |
50ff0bcd | 1724 | |
fee9ef0d | 1725 | return fCurrentVertex; |
eae1677e | 1726 | } |
50ff0bcd | 1727 | //---------------------------------------------------------------------------- |
6d8df534 | 1728 | AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray, |
97669588 | 1729 | Bool_t optUseFitter, |
1730 | Bool_t optPropagate, | |
1731 | Bool_t optUseDiamondConstraint) | |
1732 | ||
fee9ef0d | 1733 | { |
50ff0bcd | 1734 | // |
6d8df534 | 1735 | // Return vertex from array of ESD tracks |
50ff0bcd | 1736 | // |
5834e9a1 | 1737 | |
6d8df534 | 1738 | Int_t nTrks = (Int_t)trkArray->GetEntriesFast(); |
1739 | UShort_t *id = new UShort_t[nTrks]; | |
1740 | ||
1741 | AliESDtrack *esdt = 0; | |
50ff0bcd | 1742 | for(Int_t i=0; i<nTrks; i++){ |
6d8df534 | 1743 | esdt = (AliESDtrack*)trkArray->At(i); |
1744 | id[i] = (UShort_t)esdt->GetID(); | |
50ff0bcd | 1745 | } |
1746 | ||
97669588 | 1747 | VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint); |
6d8df534 | 1748 | |
4ce766eb | 1749 | delete [] id; id=NULL; |
6d8df534 | 1750 | |
1751 | return fCurrentVertex; | |
50ff0bcd | 1752 | } |
1753 | //-------------------------------------------------------------------------- |